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ABSTRACT 

A probabilistic technique for the joint estimation of background and sources with the 
aim of detecting faint and extended celestial objects is described. Bayesian probabil- 
ity theory is applied to gain insight into the coexistence of background and sources 
through a probabilistic two-component mixture model, which provides consistent un- 
certainties of background and sources. A multi-resolution analysis is used for revealing 
faint and extended objects in the frame of the Bayesian mixture model. All the re- 
vealed sources are parameterized automatically providing source position, net counts, 
morphological parameters and their errors. 

We demonstrate the capability of our method by applying it to three simulated 
datasets characterized by different background and source intensities. The results of 
employing two different prior knowledges on the source signal distribution are shown. 
The probabilistic method allows for the detection of bright and faint sources indepen- 
dently of their morphology and the kind of background. The results from our analysis 
of the three simulated datasets are compared with other source detection methods. 
Additionally, the technique is applied to i?05'^T' all-sky survey data. 
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1 INTRODUCTION 

Background estimation is an omnipresent problem for source 
detection methods in astrophysics, especially when the 
source signal is weak and difficult to discriminate against 
the background. An inaccurate estimation of the background 
may produce large errors in object photometry and the loss 
of faint and/or extended objects. 

Deep observations are commonly obtained by combin- 
ing several individual exposures to generate a final astro- 
nomical image. Often large exposure variations characterize 
these data and the background may vary significantly within 
the field. Hence, the background modelling has to incorpo- 
rate the knowledge provided by the observatory's exposure 
time without compromising the statistical properties. In ad- 
dition, instrumental structures, such as detector ribs or CCD 
gaps, produce lack of data. The missing data must be han- 
dled consistently in the background estimation to prevent 
undesired artificial effects. 

Celestial objects exhibit a large variety of morphologies and 
apparent sizes. Thus, the search for sources should not be 
driven by any predefined morphology, allowing proper esti- 
mate of their structural parameters. The instrumental point 
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spread function (PSF) is often not known exactly for the 
whole field of view. So a source detection algorithm should 
be able to operate effectively without the knowledge of the 
instrument characteristics Q. 

Conventional source detection met hods employed in 
deep imaging surveys are: SExtract or jBertin fc Arnoutd 
| l996|), slid i ng w indow technique l|Harnden et al.l Il984l : 
Gioia et all ll990D and maximum likelihood PSF fit- 



ting fHasing er et all Il994l : 
wavelet transformation (e.. 



Rosati et al, 
19981 : 



Il995l: 



Freeman et al 



Boese fc Doebereincr 2001 



Slezak. Biiaoui fc Mars 199^ 



5' 



Damiani et al.l Il997l ; IStarck fc Pierrel 
"2002!). A review of t hese techniques 
can be found in IValtchanov. Pierre fc Gast aud (2001) and 
iBecker et al] (|2007l ). 

SExtractor is one of the most widely used source detec- 
tion procedure in astronomy. It has a simple interface and 
very fast execut ion. It p roduces reliable aperture photom- 
etry catalogues l| Becker et al. 200^). It is applied in X-ray 
regime on filtered images (Valtchanov et al. 2001). 

The sliding window technique is a fast and robust source 



^ If the instrumental PSF is known precisely, then that informa- 
tion can be taken into account in a further step of source charac- 
terization within the properties of a specific observational setup. 



2 F. Guglielmetti et al. 



detection method. This technique may fail while detect- 
ing extended sources, sources near the detection limit and 
nearby sources (Valtchanov et al. 2001). This source detec- 
tion method has been refined with more elaborated tech- 



niques, such as matched filters (e.g. IVikhliiiin et alj 19951 : 



IStewartl200l ) and recently the Cash method ijStewar tl2009l l. 
The Cash method is a maximum likelihood technique. For 
source detection, the method employs a Cash likelihood- 
ratio stat istic, that i s an extended statistic for Pois- 
son data (|Cashlll979h . Both the matched filters and Cash 
methods are at least by a fact or of 1.2 mor e sensitive than 
the sliding-window technique l|Stewartll2009l ). Though, both 
methods are designed for the detection of point sources. 
The candidate sources are characterized in a further step 
using maximum likelihood PSF fitting. The maximum like- 
lihood PSF fitting procedure performs better than other 
conventional techniques for flux measurements of point-like 
sources although accurate photometry is achieved when the 
PSF model used is close to the image PSF (Valtchanov et 
al. 2001). In iPierre et"aLl (|2004l '), the maximum likelihood 
profile fit on photon images is extended taking into account 
a spherically symmetric /3-model (King profile, see refs.JKing, 
1 19621 . ICavaliere fc Fusco-Femiano llQTSl ) convolved with the 
instrumental PSF for improving the photometry of extended 
objects. 

Wavelet transform (WT) techniques improve the de- 
tection of faint and extend ed sources with respec t to other 
conventional methods (see IStarck fc Pierre! Il998l for more 
details). In fact, WT techniques are able to discriminate 
structures as a function of scale. Within larger scales, faint 
and extended sources are revealed. WTs are therefore valu- 
able tools for the detection of both point -like and extended 
sources (Valtchanov et al. 2001). Nonetheless, these tech- 
niques favor the detection of circularly symmetric sources 
(Valtchanov et al. 2001). In addition, artefacts may appear 
around the detected structur es in the reconstructed image, 
and the flux is not preserved jStarck fc Pierrelll998l ). In or- 
der to overcome these problems, some reflnement s have been 
applied to the WT techniques. IStarck fc Pierre! (11998! ). for 
instance, employ a multi-resolution support filtering to pre- 
serve the flux and the adjoint WT operator to suppress arte- 
facts which may appear aro und the objects. An advance on 
this method is presented in !Starck et al.! (!2000! ). A space- 
variant PSF is incorporated in their WT technique. Object 
by object reconstruction is performed. For point sources the 
flux measurements are close to that obtained by PSF fitting. 

The detection of faint extended and point-like sources 
is a non-trivial task for source detection methods. This task 
becomes more complicated when the detection rate and pho- 
tometric reconstruction for the sources are taken into ac- 
count (Valtchanov et al. 2001). 

A self-consistent statistical approach for background 
estimation and source detection is given by Bayesian 
probability theory (BPT), which provides a general and 
consistent frame for logical inference. The achievement of 
Bayesian techniques on signal detections in astrophysics 
has already been s hown , f or example in the w o rks o f 
iGregorv fc Loredd (|l992! ). iLoredo fc Wasserman! (!l995! ) 
and !Scargle! for instance. In modern observa- 

tional astrophysics, BPT techniqu es for image analysis 
have b een extensively applied: e.g. iHobson fc McLachlan 
l!2003! ). [Caxvalho. Rocha fc HobsonI (|2008! ). !Savage fc Oliver 



|2003), !Strong! 12003! '). For the detection of discrete ob- 
jects embedded in G a ussian noise (microwave regime), 
iHobson fc McLachlan! (|2003! ') utilizes a model-fitting 
methodology, where a parameterized form for the objects of 
interest is assumed. Markov-chain Monte Carlo (MCMC) 
techniques are used to explore the parameter space. An 
advance on this work is provided by Carvalho et al. (2008). 
For s peeding up the method of IHobson fc McLachlan! 
(|2003! ). Carvalho et al. (2008) proposes to use Gaussian 
approximation to the posterior probability density function 
(pdf) peak when performing a B ayesian model se l ection 
for source detection. The work of ISavage fc Olh^ (|2007! ) 
is developed within Gaussian statistics (infrared data). At 
each pixel position in an image, their method estimates the 
probability of the data being described by point source or 
empty sky under the assumptions that the background is 
uniform and the sources have circular shapes. The Bayesian 
information criterion is used for the selection of the two 
models. Source parameters are estimated in a second 
step employing Bayesian model selection. !Strond (|2003! ) 
developed a technique for image analysis within Poisson 
statistics. The technique is instrument specific and is 
applied to 7-ray data. The first objective of this technique 
is to reconstruct the intensity in each image pixel given 
a set of data. The Maximum Entropy method is used for 
selecting from the data an image between all the available 
ones from a multi-dimensional space. The dimension of the 
space is proportional to the number of image pixels. 

We propose a new source detection method based on 
BPT combined with the mixture-model technique. The al- 
gorithm allows one to estimate the background and its un- 
certainties and to detect celestial sources jointly. The new 
approach deals directly with the statistical nature of the 
data. Each pixel in an astronomical image is probabilisti- 
cally assessed to contain background only or with additional 
source signal. The results are given by probability distri- 
butions quantifying our state of knowledge. The developed 
Background-Source separation (BSS) method encounters: 
background estimation, source detection and characteriza- 
tion. 

The background estimation incorporates the knowledge of 
the exposure time map. The estimation of the background 
and its uncertainties is performed on the full astronomi- 
cal image employing a two-dimensional spline. The spline 
models the background rate. The spline amplitudes and the 
position of the spline supporting points provides fiexibil- 
ity in the background model. This procedure can describe 
both smoothly and highly varying backgrounds. Hence, no 
cut out of regions or employment of meshes are needed 
for the background estimation. The BSS technique does 
not need a threshold level for separating the sources from 
the background as conventional methods do. The thresh- 
old level is replaced by a measure of probability. In conven- 
tional methods, the threshold level is described in terms of 
the noise standard deviation, then translated into a prob- 
ability. The classification assigned to each pixel of an as- 
tronomical image with the BSS method allows one to de- 
tect sources without employing any predefined morpholo- 
gies. Only, for parametric characterization of the sources 
predefined source shapes are applied. The estimation of 
source parameters and their uncertainties includes the es- 
timated background into a forward model. Only the statis- 
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tics of the original data are taken into account. The BSS 
method provides simultaneously the advantages of a multi- 
resolution analysis and a multi-color detection. Objects in 
astronomical images are o rganized in hierarchical structures 
l|Starck fc Murtaghl [200^ ') . In order to quantify the multi- 
scale structur e in the data, a multi-resolution ana,lysis i s 
required (see IStarck et all I2OO0I : iKolaczvk fc DixonI I2OO0I ) . 
In the BSS approach the multi-resolution analysis is incor- 
porated in combination with the source detection and back- 
ground estimation technique with the aim to analyse statis- 
tically source structures at multiple scales. When multiband 
images are available, the information contained in each im- 
age can be statistically combi ned in order to extend the 
detection limit of the data (see Szalav^ Connollv fc Szokolvl 



Bernardo fc Smith! Il994l: ICelman et all Il995l: [Shdal Il996l: 



ll999l : lMurtagh. Rafterv fc Starckll200a )" 

The capabilities of this method are best shown with the de- 
tections of faint sources independently of their shape and 
with the detections of sources embedded in a highly vary- 
ing background. The technique for the joint estimation of 
background and sources in digital images is applicable to 
digital images collected by a wide variety of sensors at any 
wavelength. The X-ray environment is particularly suitable 
to our Bayesian approach for different reasons. First, X- 
ray astronomy is characterized by low photon counts even 
for relatively large exposures and the observational data are 
unique, i.e. the experiment is rarely reproduced. Second, the 
X-ia.y astronomical images provided by new generation in- 
struments are usually a combination of several individual 
CCD imaging pointings. Last, there are few source detec- 
tion algorithms developed so far for an automated search of 
faint and extended sources. 

In this paper, our aim is to describe the developed BSS 
technique. The outline of this paper is as follows. In Section 
[2I we briefly review the basic aspects of BPT. In Section O 
we introduce the background estimation and source detec- 
tion technique. In Section O the BSS algorithm is extended 
in order to obtain an automated algorithm for source charac- 
terization. In Section [5] the issue of false positives in source 
detection is addressed. In Section |S1 the BSS method is ap- 
plied to simulated data. We show results for two different 
choices of prior pdf for the source signal. In Section!?] our re- 
sults on the three simulated datasets are compared with th e 
outcome from wavdetect algorithm (jPreeman et al.ll2002h . 
In Section [8] we test the BSS method on astronomical im- 
ages coming from ROSAT all-sky survey (RASS) data. Our 
conclusions on the BSS technique are provided in Section |9l 
In Section [TOl a list of acronyms used throughout the paper 
is reported. 

The flexibility of our Bayesian technique allows the in- 
vestigation of data coming from past and modern instru- 
ments without changes in our algorithm. An extensive ap- 
pUcation of the BSS technique to real data is addressed in a 
forthcoming paper (Guglielmetti et al., in preparation). In 
the latter, the BSS algorithm is used for the analysis of a 
specific scientific problem, i.e. the search of galaxy clusters 
in the CDF-S data. 



2 BAYESIAN PROBABILITY THEORY 

In order to analyse the heterogeneous data present in as- 
tronomical images, we employ BPT (see Ijeffrevsl Il96ll : 



Javnesll2003l : lDosJl200alO'Hagan fc Forsteij|2004 ICregorvl 
20051). BPT gives a calculus how to resolve an inference 
problem based on uncertain information. The outcome of 
the BPT analysis is the pdf of the quantity of interest, 
which encodes the knowledge to be drawn from the infor- 
mation available (a posteriori). The posterior pdf comprises 
the complete information which can be inferred from the 
data and the statistical model, and supplementary informa- 
tion such as first-principle physics knowledge. 

This statistical approach is based on comparisons 
among alternative hypotheses (or models) using the sin- 
gle observed data set. Each information entering the models 
that describe the data set is combined systematically. The 
combination includes all data sets coming from different di- 
agnostics, physical model parameters and measurement nui- 
sance parameters. Each data set and parameters entering 
the models are subject to uncertainties which have to be 
estimated and encoded in probability distributions. Within 
BPT the so-called statistical and systematic uncertainties 
are not distinguished. Both kinds of uncertainties are treated 
as a lack of knowledge. 

Bayes' theorem states: 



PiHi\D,a,I) = 



P{D\Huo, I)Pmi) 

Pirn 



(1) 



where the number of alternative hypotheses to be compared 
is larger than or equal to 2. Bayes' theorem is a consequence 
of the sum and product rules of probability theory. The ver- 
tical bars in ^ denote conditionality property, based on 
either empirical or theoretical information. Equation ^ re- 
lates the posterior pdf P{Hi\D, a, I) to known quantities, 
namely, the likelihood pdf P{D\Hi, a, I) and the prior pdf 
P{Hi\I). P{D\I) is the evidence of the data which consti- 
tutes the normalization and will not affect the conclusions 
within the context of a given model. The posterior pdf is the 
quantity to be inferred. It depends on the full data set D, 
on the errors a entering the experiment and on all relevant 
information concerning the nature of the physical situation 
and knowledge of the experiment (7). The likelihood pdf 
represents the probability of finding the data D for given 
quantities of interest, uncertainties a and additional infor- 
mation /. It reveals the error statistics of the experiment. 
The prior pdf is independent from the experimental errors. 
It represents physical constraints or additional information 
from other diagnostics. The terms 'posterior' and 'prior' 
have a logical rather than temporal meaning. The poste- 
rior and prior pdfs can be regarded as the knowledge 'with' 
and 'without' the new data taken into account, respectively. 

In order to arrive at the pdf of any quantity x in the 
model, marginalization of the multi-dimensional pdf can be 
regarded as a projection of the complete pdf on to that quan- 
tity. Marginalization is performed by integration over the 
quantity y one wants to get rid of: 



(2) 



P{x)^ / P{x,y)dy= / P{x\y)P{y)dy 



Marginalization of a quantity y thus takes into account the 
uncertainty of y which is quantified by the pdf P(y). The un- 
certainty of y propagates into the pdf P(x). Marginalization 
provides a way to eliminate variables which are necessary to 
formulate the likelihood but otherwise uninteresting. 
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Another important property of BPT is the capability 
of modelling t he data by mixture di s tributions in the pa- 
rameter space (|Everitt fc Handlll98ll : lNeallll992l ). Mixture 
distributions are an appropriate tool for modelling processes 
whose output is thought to be generated by several different 
underlying mechanisms, or to come from several different 
populations. Our aim is to identify and characterize these 
underlying 'latent classes'. Therefore, we follow the stan- 
dard Bayesian approach to this problem, which is to define 
a prior distribution over the parameter space of the mix- 
ture model and combine this with the observed data to give 
a posterior distribution over the parameter space. Essential 
for model comparison or object classification is the marginal 
likelihood ( evidence, prior predicted value) . Marginalization 
(integration) of the likelihood over parameter space provides 
a measure for the credibility of a model for given data. Ra- 
tios of marginal likelihoods (B ayes factors) are frequ ently 
used for comparing two models ()Kass fc Raftervl[l995h . De- 
tails of mixture modelling in the framework of BPT can be 
found in von der Linden et al. (199 7, 1999) and Fis c her et 
al. (2000, 2001, 2002). In particular. [Fischer fc d"os3 \2QQ± 
have demonstrated the capability of the Bayesian mixture 
model technique even with an unknown number of compo- 
nents for background separation from a measured spectrum. 
The present approach follows these previous works. 



3 THE JOINT ESTIMATION OF 

BACKGROUND AND SOURCES WITH BPT 

The aim of the BSS method is the joint estimation of back- 
ground and sources in two-dimensional image data. We can 
identify two basic steps of our algorithm: (A) background 
estimation and source detection, (B) calculation of source 
probability maps in a multi-resolution analysis. 

The input information of the developed algorithm is the 
experimental data, i.e. the detected photon counts, and the 
observatory's exposure time. 

The background rate is assumed to be smooth, e.g. spa- 
tially slowly varying compared to source dimensions. To al- 
low for smoothness the background rate is modelled with a 
two-dimensional thin-plate spline (TPS), ('Section l3.2p . The 
number and the positions of the pivots, i.e. the spline's sup- 
porting points, decide what data structures are assigned to 
be background. All structures which can not be described 
by the background model will be assigned to be sources. The 
number of pivots required to model the background depends 
on the characteristics of the background itself. Though the 
minimum number of pivots is four (since the minimum ex- 
pansion of the selected spline has four terms), their number 
increases with increasing background variation. 

The coexistence of background and sources is described 
with a probabilistic two-component mixture model (Section 
I3.1|l where one component describes background contribu- 
tion only and the other component describes background 
plus source contributions. Each pixel is characterized by the 
probability of belonging to one of the two mixture compo- 
nents. For the background estimation the photons contained 
in all pixels are considered including those containing addi- 
tional source contributions. No data censoring by cutting 
out source areas is employed. 

For background estimation the source intensity is con- 



sidered to be a nuisance parameter. According to the rules 
of BPT, the source signal distribution is described prob- 
abilistically in terms of a prior pdf. The prior pdf of the 
source signal is an approximation to the true distribution of 
the source signal on the field. We studied two prior pdfs of 
the source signal: the exponential and the inverse-Gamma 
function. 

The background and its uncertainties (Section 13. 3|) are 
estimated from its posterior pdf. Therefore, for each pixel 
of an astronomical image an estimate of its background and 
its uncertainties are provided. 

Moreover, the Bayesian approach introduces hyper- 
parameters, that are fundamental for the estimation of 
the posterior pdfs for the background and source intensi- 
ties. Specifically, in Section 13.41 we show that the hyper- 
parameters are estimated exclusively from the data. 

The source probability is evaluated with the mixture 
model technique for pixels and pixel cell^ in order to en- 
hance the detection of faint and extended sources in a multi- 
resolution analysis. Pixels and pixel cells are treated iden- 
tically within the Bayesian formalism. For the correlation 
of neighbouring pixels, the following methods have been in- 
vestigated: box filter with a square, box filter with a circle, 
Gaussian weighting filter (see Section [3. II for more details). 

The BSS technique is morphology free, i.e. there are no 
restrictions on the object size and shape for being detected. 
An analysed digital astronomical image is converted into the 
following: 

I) the background rate image, or 'TPS map', is an array 
specifying the estimated background rate at each image pixel 
for a given observation. The effects of exposure variations are 
consistently included in the spline model; 
II) the background image, or 'background map', is an array 
specifying the estimated background amplitude at each im- 
age pixel for a given observation. It is provided by the TPS 
map multiplied with the telescope's exposure; 
III) the source probability images, or 'source probability 
maps' (SPMs), display the probability that source counts 
are present in pixels and pixel cells for a given observation 
in a multi-resolution analysis. 

Movies are produced with the SPMs obtained at different 
resolutions. The moving images allow one to discern inter- 
actively the presence of faint extended sources in digital 
astronomical images. The size of faint extended sources is 
correlated with the scale of the resolution, used for their 
detection. SPMs coming from other energy bands can be 
combined statistically to produce conclusive SPMs at dif- 
ferent resolutions with the advantage to provide likelihoods 
for the detected sources from the combined energy bands 
(Section EXg. 

3.1 Two— component mixture model 

The general idea of the described Bayesian model is that a 
digital astronomical image consists of a smooth background 
with additive source signal, which can be characterized by 

^ We define pixels as the image finest resolution Hmited by instru- 
mental design, while we define pixel cells as a group of correlated 
neighbouring pixels. 
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any shape, size and brightness. The background signal is 
the diffuse cosmic emission added to the instrumental noise. 
The source signal is the response of the imaging system to 
a celestial object. A surface 6(x) describes the background 
under the source signal, where x = (x, y) corresponds to the 
position on the grid in the image. 

Therefore, given the observed data set D — {dij} G No, 
where dij is photon counts in pixel (or pixel cell) {ij}, two 
complementary hypotheses arise: 



dii = h + Eij 

d 



bij + Sij + eij. 



Hypothesis Bij specifies that the data consists only of 
background counts bij spoiled with noise tij, i.e. the (statis- 
tical) uncertainty associated with the measurement process. 
Hypothesis Bij specifies the case where additional source 
intensity Sij contributes to the background. Additional as- 
sumptions are that no negative values for source and back- 
ground amplitudes are allowed and that the background is 
smoother than the source signal. 

The smoothness of the background is achieved by modelling 
the background count rate with a bivariate TPS where the 
supporting points are chosen sparsely to ensure that sources 
can not be fitted. The spline fits the background component 
whereas count enhancements classify pixels and pixel cells 
with source contributions. 

In the following, pixel cells are subsumed by pixels. 
Pixel cells are collections of pixels where dij is the total pho- 
ton count in the cell {ij} - The photon counts of neighbouring 
pixels are added up and the formed pixel cell is treated as 
a pixel. In principle, any cell shape can be chosen. In prac- 
tice, two methods have been developed when pixels have 
weight of one within the cell (box filtering with cell shape 
of a square or of a circle) and one method when pixels have 
different weights within the cell (Gaussian weighting). The 
box filter with cells of squared shape consists of taking in- 
formation of neighbouring pixels within a box. The cell size 
is the box size. The box filter with cells of circular shape 
considers pixels with a weight of one if inside the cell size, 
otherwise zero. Pixels have a weight of one when the cell size 
touches them at least at the centre. This method allows the 
pixel cells to have varying shapes. The Gaussian weighting 
method provides Gaussian weights around a centre: weights 
are given at decreasing values according to the distance from 
the centre. 

As expressed in ((T}, we are interested in estimating the 
probabilities of the hypotheses: Bij and _Bij. 
In this paper we address the problem when the photon 
counts are small and Poisson statistics have to be used. The 
likelihood probabilities for the two hypotheses within Pois- 
son statistics are: 



by , 

p(dijlBij,&ij) = -^e-''^ 

p(d,|Bi,,6i,,.ij) = i^li±£iil^ie-(^+-^). 

dij! 



(3) 
(4) 



This technique is easily adaptable to other likelihoods, 
for instance emplo ying Gaussian statistics as given in 
iFischer et all (j200lh . 

The prior pdfs for the two complementary hypotheses are 
chosen to be p(-Bij) = /3 and p(-Bij) = 1-/3, independent of 
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Figure 1. Representation of selected prior pdfs of the source sig- 
nal versus the photon counts. The exponential prior pdf is drawn 
at two A values (dashed lines). The inverse-Gamma function 
prior pdf is plotted at two a values (continuous lines). On the 
ordinate, 7 indicates A or a. 



i and j. Specifically, the parameter (3 is the prior probability 
that a pixel contains only background. 

Since it is not known if a certain pixel (or pixel cell) 
contains purely background or additional source signal, the 
likelihood for the mixture model is employed. The likelihood 
for the mixture model effectively combines the probability 
distributions for the two hypotheses, Bij and Bij-. 

p(D|b,s) = 

= n^^ll^"'^ + (1 - P) ^''^^f' e-''''+^^'^, (5) 
ij 

where b — {bij}, s = {sij} and {ij} corresponds to the pixels 
of the complete field. 

The probability of having source contribution in pix- 
els and p ixel cells is according to Bayes' theorem (see 
iFischer et al. 2001): 



P{Bij 
(1^ 



dij , bij , Sij ) — 



-f3)-p{dij 



Bij, bij, Sij) 



/3-p(dij 



j) + (l-/3)-p(dij I Bij, bij, Sij) 



(6) 



This equation enables the data contained in an astronomical 
image to be classified in two groups: with and without source 
signal contribution. 

Equation (|6} is also used in the multi-resolution analy- 
sis. The SPM with the largest resolution is characterized by 
the probability of uncorrelated pixels. At decreasing resolu- 
tions a correlation length is defined. Starting from a value of 
0.5, the correlation length increases in steps of 0.5 pixel for 
decreasing resolution. The SPMs at decreasing resolutions 
are, therefore, characterized by the information provided 
by background and photon counts in pixel cells. Specifi- 
cally, photon counts and background counts are given by 
a weighted integration over pixel cells. The integrated pho- 
ton and background counts enter the likelihood for the mix- 
ture model. Then, the source probability is estimated for 
each image pixel in the multi-resolution analysis. The multi- 
resolution algorithm preserves Poisson statistics. 
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3.1.1 Source signal as a nuisance parameter 

Following Fischer et al. (2000, 2001), the source signal in @ 
is a nuisance parameter, which is removed by integrating it 
out (marginalization rule, eq. ([2}): 



p(dij|Bij,feij,Sij)-p(sij|Bij,6ij)dsij. 



A nuisance parameter is a parameter which is important for 
the model (describing the data), but it is not of interest at 
the moment. Following BPT, a prior pdf of the source signal 
has to be chosen. The final result depends crucially on the 
prior pdf set on the source signal. In fact, in addition to the 
choice of the TPS pivots, the prior pdf of the source signal 
provides a description of what is source. All that is not de- 
scribed as a source is identified as background and vice versa. 
Two approaches are presented: the first method accounts 
for the knowledge of the mean value of the source intensity 
over the complete field {exponential prior), the second ap- 
proach interprets the source signal distribution according to 
a power-law (inverse-Gamma function prior). 

3.1.1.1 Exponential prior Following the works of Fis- 
cher et al. (2000, 2001, 2002), we choose a prior pdf of the 
source signal that is as weakly informative as possible. The 
idea follows a practical argument on the difficulty of provid- 
ing sensible information. We describe the prior pdf on the 
source intensity by an exponential distribution. 



p(sij I A) 



A 



(7) 



This is the most uninformative Maximum Entropy distribu- 
tion for known mean value of the source intensity A over the 
complete field. 

In Fig. [1] equation ^ is drawn for two values of the mean 
source intensity: A = 1 count and A = 10 counts. No bright 
sources are expected to appear in fields with A ~ 1. In the 
case of values of A ^ 1, bright and faint sources are present 
in these fields. 

The marginal Poisson likelihood for the hypothesis Bij has 
the form: 



p(dij|Bij,6ij, A) = 



hi 



A(l-K i)(''.j+i) 

^ r[(dij + i),&ij(i + i)] 



(8) 



r(dij + 1) 

where r[o, a;] = e~^t°'~^dt (a > 0) is the incomplete- 
Gamma function and T[a] = r[a,0] is the Gamma function 
(see refs. Fischer et al. 2000, 2001). 

The behaviour of the Poisson and the marginal Poisson 
distributions is depicted with a parameter study. For the 
parameter study three background amplitudes b are chosen: 
0.1, 1 and 10 counts. In Fig. [2]the Poisson distribution ([3]) 
and the marginal Poisson distribution ([SJ are drawn on a 
logarithmic scale. These plots are indicated with (a), (6) and 
(c), respectively. The parameter A, which describes the mean 
intensity in a field, has values of: 1, 10 and 100 counts. The 
selected values for the background and for the parameter A 
are chosen such that the large variety of cases one encoun- 
ters analysing digital astronomical images are covered. For 
instance, 6 = 10 counts and A = 1 count (plot (c)) corre- 
sponds to a field when the source signal is much smaller than 



the background signal. On the other side, 6 = 0.1 counts and 
A = 100 counts (plot (a)) corresponds to a field character- 
ized by bright sources and small background amplitude. 

The Poisson distribution is larger than the marginal 
Poisson distribution for photon counts lower than or equal 
to the background intensity. Hence, hypothesis Bij is more 
likely than the complementary hypothesis Bij. This situa- 
tion changes when the photon counts are larger than the 
background amplitude. 

The decay length of the marginal Poisson distribution 
is determined by the expected source intensity A. The prob- 
ability to detect pixels satisfying hypothesis B^j is sensitive 
to the decay length of the marginal Poisson distribution and 
to the background amplitude, that is recognizable in the dis- 
tance between the Poisson and the marginal Poisson distri- 
butions. Hence, the BSS method allows probabilities to be 
sensitive to the parameters characterizing the analysed dig- 
ital astronomical image. 

Let's consider plot (b) in Fig. [2] for photon counts in the 
range (0 — 10). The background amplitude has a value of 1 
count. If the expected mean source intensity on the complete 
field has a value of 1 count, i.e. A = 1 count, 3 photon counts 
or more in a pixel are classified as a source. The probability 
of detecting a source increases with increasing counts in a 
pixel. This is due to the increasing distance of the marginal 
Poisson likelihood from the Poisson likelihood. If an analyst 
allows for many bright sources distributed in the field, then 
the relative number of faint sources is reduced. In fact, when 
a mean source signal 100 times larger than the background 
is expected, then 5 photon counts or more in a pixel are 
needed to identify the event due to the presence of a source. 
When A — 100 counts, 5 photon counts in a pixel reveal a 
source probability lower than the one obtained when A = 1 
count. This situation changes for 7 or more photon counts 
in a pixel. 

3.1.1.2 Inverse— Gamma function prior Alterna- 
tively to the exponential prior, a power-law distribution can 
be chosen to describe the prior knowledge on the source sig- 
nal distribution. 

We are tempted to claim that any physical situa- 
tion contains sensible information for a proper (normal- 
izable) prior. We choose a prior pdf that is inspired by 
the cumulative signal number counts distribution used of- 
ten in astrophysics for describing the integral fiux distribu- 
tion of cosmological sources, i.e. a power-law function (see 
refs. iRosati. Borgani fc NormarJ |2002| : iBrandt fc Hasingeij 
I2OO5I and references therein) . The power law can not be em- 
ployed as a prior pdf, because the power law is not normal- 
ized. We use instead a normalized inverse-Gamma function. 
It behaves at large counts as the power law, because it is 
described by a power law with an exponential cutoff. 

The prior pdf of the source signal, described by an 
inverse-Gamma function, is: 



p(sij 1 a, a) 



r(a-l)' 



a > 1; a > 0, 



(9) 



with slope a and cutoff parameter a. When a has a small 
positive value, the inverse-Gamma function is dominated by 
a power-law behaviour. 

The parameter a gives rise to a cutoff of faint sources. 
This parameter has important implications in the estima- 
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Figure 2. Likelihood pdfs versus photon counts. The Poisson distribution is indicated with PD. The marginal Poisson distribution is 
indicated with MPD. Plots (a), (fe) and (c) are obtained using the exponential prior pdf. Plots (d), (e) and (/) are created using the 
inverse-Gamma function prior pdf. 



tion of the background. If a is smaller than the background 
amplitude, the BSS algorithm detects sources of the order 
of the background. If a is larger than the background ampli- 
tude, the BSS algorithm assigns faint sources with intensities 
lower than a to be background only. 

In Fig. [T] equation (|9]) is drawn for two values of the 
parameter a. For this example the cutoff parameter a has a 
value of ~ 0.1 count. The distributions peak around a. The 
decay of each distribution depends on the value of a. When 
a is large, i.e. for values > 2, the distribution drops quickly 
to zero. Instead the distribution drops slowly to zero, when 
a approaches 1. Hence, small values of a indicate bright 
sources distributed on the field. 

The marginal Poisson likelihood for the hypothesis Bij 
is now described by: 



p{dii\Bij,bij,a,a) = 



2 -i,, 
e 'J 



r(a-i) 



■r(fc + i)r(dij -fc + 1) 



(10) 



where a > 1, a > and Kk-a+ii'2\'^) is the Modified Bessel 
function of order k — a + 1. 

In order to avoid numerical problems with the Bessel func- 
tion, the following upward recurrence relation was derived: 



{k-l)k z-k 



where K^{z) ■- K^{z)/T{k + 1) and v ^ k - a + 1. K^{z) 
is the Bessel function of imaginary argument and it has the 
property: K-i,{z) = K^{z). 

Fig- [21 {d)-{f), shows the corresponding parameter 
study for the inverse-Gamma function prior. The parameter 




5 10 15 

photon counts 



20 



Figure 3. Likelihood pdfs for the mixture model using the ex- 
ponential prior (dashed line) and the inverse-Gamma function 
prior (continuous line). The Poisson and the marginal Poisson 
likelihood pdfs are weighted with (3 and (1 — /3), respectively. 



a is chosen to be 1.3, 2.0 and 3.0 and the cutoff parameter a 
~ 0.1 counts. The decay length of the marginal Poisson dis- 
tributions (|10p are now given by the value of a. The decay 
length decreases with increasing a values. 



3.1.2 The likelihood for the mixture model 

The marginal Poisson likelihood pdf will be indicated with 
p(dij|-Bij, 7), where 7 indicates A or a. In the case of the 
inverse-Gamma function prior pdf, the cutoff parameter a 
does not appear since the value of this parameter is chosen 
such that the inverse-Gamma function is dominated by a 
power-law behaviour. 

The likelihood for the mixture model, as written in ([5]), 
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now reads: 



p(D i 6,7,/3) = \[[l3-p{d,^ i Bij,feij) + (1 - /3) 



xp(dij I Bij,6ij,7)]; D = {dij}, & = {&ij}. 



(11) 



In Fig.|3]we show the effect of the hkelihood for the mixture 
model on the data (semi-log plot). The likelihood pdf for 
the mixture model is drawn for each prior pdf of the source 
signal employing the background value b and the prior pdf (3. 
The value chosen for the parameter (5 indicates that 99 per 
cent of the pixels distributed in the astronomical image are 
containing only background. The likelihood pdfs are com- 
posed by a central peak plus a long tail. The central peak 
is primarily due to the presence of the Poisson distribution. 
The long tail is caused by the marginal Poisson distribution. 
The presence of the long tail is essential in order to reduce 
the influence of the sou rce signal for background estimation 
(see iFischer et al]|2000l for more details). 




3.2 Thin— plate spline 

The TPS has been selected for modelling the structures aris- 
ing in the background rate of a digital astronomical image. 
It is a type of radial basis function. 

The TPS is indicated with t(x), where x = {x,y) corre- 
sponds to the position on the grid in the detector field. The 
shape of the interpolating TPS surface fulfills a minimum 
curvature condition on infinite support. 

More specifically, the TPS is a weighted sum of trans- 
lations of radially symmetric basi s functions augmented by 
a linear term (see lMeinguet|[l979l and IWahba.2006,' ). of the 
form 



t(x) 



i5(x)+^Ai/(x- 



xi), X e R 



-E(x) = Co + cix + C2y is the added plane. is the number 
of support points (pivots). The weight is characterized by 
Al. /(x — xi) is a basis function, a function of real values 
depending on the distance between the grid points x and 
the support points xi, such that |x — xi| > 0. 
Given the pivots xi and the amplitude Zi = ^(xi), the TPS 
satisfies the interpolation conditions: 



t(xi) = i5(xi) 



and minimizes 



E 



Ai/(xi - xi) 



(./xx + 2/xy + fyy)dxdy. 



\\t\\'^ is a measure of energy in the second derivatives of t. In 
other words, given a set of data points, a weighted combi- 
nation of TPSs centered about each pivot gives the interpo- 
lation function that passes through the pivots exactly while 
minimizing the so-called 'bending energy'. The TPS satisfies 
the Euler-Lagrange equation and its solution has the form: 
/(x — xi) ~ r^ln(r^), where = (a; — xi)^ + {y — yi)^- This 
is a smooth function of two variables defined via Euclidean 
space distance. In Fig. [3] an example of TPS with one pivot 
is pictured. 

In order to fit the TPS to the data, it is necessary to solve 



Figure 4. Example of thin-plate spline, /(x — xj) = r^ln(r'^) is 
drawn with one support point centered at tlie Cartesian origin. 
The projection of /(x — xi) on the {x,y) plane is placed below 
and its contour plot on top of the surface image. 



for the weights and the plane's coefficients so that it is pos- 
sible to interpolate the local TPS's amplitude: 

tii = tij(iVr,{xi,2l,Z = l,...,iV,}) 

which is the background rate, foij will indicate the local back- 
ground amplitude, i.e. the multiplication of tij and the local 
value of the telescope's exposure time (r): 

''ij ~ ^ '''ij- 

The TPS interpolant is defined by the coefficients, ci of the 
plane -E(x) and the weights Ai of the basis functions. The 
solution for the TPS has been evaluated on an infinite sup- 
port, since no solutions exist on a finite support, where the 
requirements for this function to be fulfilled are: 

(i) t(x) is two times continuously differentiable, 

(ii) t(x) takes a particular value Zi at the point xi, 

(iii) I[f{x,y)] is finite. 

Given the interpolation values z — {zi, . . . , zn^), the weights 
Al and ci are searched so that the TPS satisfies: 



t(xi) =21, 1= 1, 



,iVr 



and in order to have a converging integral, the following 
conditions need to be satisfied: 

e;JiAi =0 

EiJi^ij/i =0- 

This means that we have (A'r — 3) conditions on t(xi) — z\. 
The coefficients of the TPS, Ai, and the plane, a, can be 
found by solving the linear system, that may be written in 
matrix form as: 



F Q 



where the matrix components are: 
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Q 



/ I XI yi \ 

1 X2 y2 



V 1 



XN, VN, 





A 



/(Xi -Xj) 
{zi, . . . , 2jvJ^ 

(0,0,0)^ 

(Ai, . . . , Ajv,)^ 

(co,Cl,C2)^. 



After having solved (A,c)^, the TPS can be evaluated 
at any point. 

The pivots can be equally spaced or can be located 
in structure s arising in the astr ono mical imaKe. Following 
the wo rks of lFischer et~al] (|2000l ) and I von Toussaint fc Goril 
the present work can be extended employing adap- 
tive splines, i.e. allowing the number of pivots and their po- 
sitions to vary in accordance with the requirements of the 
data. 



3.3 Estimation of the background and its 
uncertainties 

The posterior pdf of the background is, according to Bayes' 
theorem, proportional to the product of the mixture like- 
lihood pdf, eq. and the prior pdf p{b), that is cho- 
sen constant for positive values of b and null elsewhere. Its 
maximum with respect to b, b* , gives an estimate of the 
background map which consists of the TPS combined with 
the observatory's exposure map. The estimation of the back- 
ground considers all pixels, i.e. on the complete field, because 
we tackle the source signal implicitly as outlier. 
The posterior pdf for b is given by: 



p{b\D)= p(b\D,z)p{z\D)d 



= J S(b-b(z)) d z. 

This integral is complicated due to the presence of the delta 
function. This is, however, of minor importance since our in- 
terest focuses on the expectation values of some functionals 
of b, say g{b). Therefore: 



E{g(b) I D) 



g{b)p(b I D)db - 



= I gibiz)f-i£l^d-'z^ 

Jg{b{z))p{D\z)piz)d'''z 
J p{D I z)p{z)d^rz 



(12) 



Assuming the maximum of p{D j z)p{z) is well determined 
we can apply the Laplace approximation: 



p{D I z)p(z)~p{D I z*)p(2*)exp{--Az^//Az}, 



(13) 



that means we approximate the integral function by a Gaus- 
sian at the maximum z* and we compute the volume un- 
der that Gaussian. The covariance of the fitted Gaussian 
is determined by the Hessian matrix, as given by (|13|) . 



where Az — 



and Hij : — 



a^lnlp(D\z)p(z)] 



is ele- 



2"'* order Taylor expansion of the posterior pdf around the 
optimized pivots amplitu de's values. For more details see 
lO'Hagan fc Forstcr' (120041 '). 
Then equation (|12l) becomes: 



E{g{b) 1 D) 



J (7(6(2))exp{-iAz^J/Az}d^ 



(27r)T^(deti/)-i 
Therefore, the posterior mean of b is: 



Eihi I D) 



VdetH 



:T/jZ* =<6ij >, 



6ij(z)exp{--!- Az^i/Az}d^'\z 



where bij{z) — Tfj ■ z, and the variance is: 



i5(A6ijA&,k I D) = I ihiiz)- < hi » 

(27r)^ 

x(6ik(2:)- < feik >)exp{-iAz^'j/Az}d^'2: 

= tT^H-^Ti^. (14) 

The la error on the estimated background function is there- 
fore calculated by the square root of equation (|14[) . 



3.4 Determining the hyper— parameters 

The two hyper-parameters 7 and /3 have so far been assumed 
to be fixed. However, these parameters can be appropriately 
estimated from the data. 

Within the framework of BPT the hyper-parameters 
(nuisance parameters) 7 and /? have to be marginalized. 
Alternatively, and not quite rigorous in the Bayesian sense, 
the hyper-parameters can be estimated from the marginal 
posterior pdf, where the background and source parameters 
are integrated out, 

maxp(/3,7 | D) ^ /3*,7*. 

Hence, the estimate of the hyper-parameters is the maxi- 
mum of their joint posterior. 

The basic idea is to use BPT to determine the hyper- 
parameters explicitly, i.e. from the data. This requires the 
posterior pdf of 7 and (3. Bayes' theorem gives: 

p{D\'y,l3)p{'y)p{P) 



p{l,P\D) = 



p{D) 



(15) 



The prior pdfs of the hyper-parameters are independent 
because the hyper-parameters are logically independent. 
These prior pdfs are chosen uninformative, because of our 
lack of knowledge on these parameters. The prior pdf for 
/3 is chosen to be constant in [0, 1]. Since 7 is a scale pa- 
rameter, the appropriate prior distribution is Jeffrey's prior: 
p(7) ~ 1/ 7. Equation (|15|l can be written as follows: 

p(7,/3 I D)(xp{D I 7,/3)p(7) = 



= P('y) / P{D \ z,"i,(3)p{z)dz 



ment {ij} of the Hessian matrix. This approximation is the 



(16) 



Assuming the maximum of p{D \ z,"f,P)p(z) is well deter- 
mined, we can apply the Laplace approximation 

p{D I z,j,l3)piz)~p{D \ z*rf,P)p{z*)exp{~^Az'^HAz} 

where Az = z — z* , z* corresponds to the maximum value 
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of the integrand in (HI, and - _ a^in\p(m.r,M-)] jg 
element {y} of the Hessian matrix. Considering dim(z) = 
A'r, where iVr is the number of support points, equation p6|) 
can be written as follows: 

(deti^)2 

p{z*) is chosen to be constant. The last term corresponds 
to the volume of the posterior pdf of 7 and {3 for each 7, j3 
estimates. 



3.5 Probability of hypothesis B 

In principle, the probability of detecting source signal in ad- 
dition to the background should be derived marginalizing 
over the background coef ficients and the hyper-param eters. 
Following the work s of Ivon der Linden et al.1 l|l999l ) and 
iFischer et all (|2000l ). the multidimensional integral, arising 
from the marginalization, can be approximated at the op- 
timal values found of the background coefficients and the 
hyper-parameters. 

Equation ^ is approximated with: 

P(Bij I dij) ~ \u..B-b') = ^— ' (17) 

^ PV°iJ l-°ij •"ij) 

T--I3 ' p(dij|"Bij,6*.,7) 

where b* = {h*j} is the estimated background amplitude, as 
explained in Section FS.SI SPMs are estimated employing this 
formula. 

The BSS method does not incorporate the shape of the 
PSF. When the PSF FWHM is smaller than the image pixel 
size, then one pixel contains all the photons coming from a 
point-like source. Otherwise point-like sources are detected 
on pixel cells as large as the PSF FWHM. Extended objects 
are detected in pixel cells large enough that the source size 
is completely included. The pixel cell must be larger than 
the PSF FWHM and it can exhibit any shape. 

Equation (|17p shows that the source probability strictly 
depends on the ratio between the Poisson likelihood, 
p{dij I -Bij,&ij), and the marginal Poisson likelihood, p(dij | 
Bijybij,'^) (Bayes factor). Bayes factors offer a way of eval- 
uating evidence in favor of competing models. 

Fig. [5] shows the effect of the mixture model technique 
on the probability of having source contribution in pixels 
and pixel cells for the exponential and the inverse-Gamma 
function prior pdfs. For the parametric studies, the param- 
eter /3 is chosen to be 0.5. This noncommittal value of (3 
arises if each pixel (or pixel cell) is equally likely to contain 
source signal contribution or background only. Psourco de- 
pends on the likelihood for the mixture model, which is the 
weighted sum of the Poisson distribution and the marginal 
Poisson distribution. For photon counts of about the mean 
background intensity, Psource is small. Psource increases with 
increasing photon counts due to the presence of the long tail 
in the marginal likelihood. This allows efficient separation of 
the source signal from the background. 

In panels (a)-(c), the distribution function of Psourco, the 
Poisson pdf (PD) and the marginal Poisson pdf (MPD) 
are drawn using the exponential prior (see also Figs. [T]and 
[2|. In the case of fields with bright sources (A > 10 times 
the background intensity), Paourco is nearly zero for photon 



Table 1. Interpretation of source probability values. 



<50% 


50% 


90% 


99% 


99.9% 


Psourco backg. 


indif. 


weak 


strong 


very 


only 








strong 



counts less than or equal to the mean background intensity. 
Psourco increases rapidly with increasing photon counts. In 
the case of fields where the mean source intensity is similar 
to the mean background intensity, pixels containing photon 
counts close to the mean background intensity have proba- 
bilities about 50 per cent. In these cases, Psourco increases 
slowly with increasing photon counts, because the two Pois- 
son distributions are similar. In the case of fields dominated 
by large background signal (A < mean background ampli- 
tude), Psourco increases very slowly with increasing photon 
counts. In this case, the decay of the marginal Poisson distri- 
bution follows closely the decay of the Poisson distribution 
(e.g. for & = 10 counts and A = 1 count). 
In Fig.Opanels (d)-(/), the distribution functions are shown 
using the inverse-Gamma function prior (see also Figs.[l]and 
[2|. The steepness of the slope depends on the parameter a 
(Fig. The source probability curve increases faster at 
decreasing a values, because small a values indicate bright 
sources distributed in the field. In panels (e) and (/), Psourco 
provides values close to 50 per cent at low numbers of pho- 
ton counts. This effect addresses the cutoff parameter a. In 
fact, in these plots the cutoff parameter a is smaller than 
the background amplitude. The situation is different in the 
simulations with small background amplitude (panel (rf)), 
where the source probabilities decrease below 50 per cent at 
low numbers of photon counts. In these simulations the cut- 
off parameter a is chosen larger than the mean background 
amplitude. Faint sources with intensities lower than a are 
described to be background. 

The interpretations of the probability of having source 
contributions in pixels and pixel cells are shown in Table 
[T] Source probabilities <50 per cent indicates the detection 
of background only. Psourco is indifferent at values of 50 per 
cent. In both cases, Psourco might contain sources but they 
can not be distinguished from the background due to statis- 
tical fluctuations. Source probability values ^ 50 per cent 
indicate source detections. False sources due to statistical 
fluctuations may occur especially for values <99 per cent. 
For more deta ils about the interpretations pr ovide d in Ta- 
ble [1] see iJeftexa (il961i ) and iKass fc Rafter^ (|l995l ). 

3. 5. 1 Statistical combination of SPMs at different energy 
bands 

An astronomical image is usually given in various energy 
bands. SPMs, obtained with (|17p . acquired at different en- 
ergy bands can be combined statistically. The probability of 
detecting source signal in addition to the background for the 
combined energy bands {fc} is: 

n 

p(Bij I dij)comb = 1 - " I ^'J)'')' (IS) 

k=l 

where n corresponds to the effective energy band. 
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Figure 5. Source probability (Psourcc)i Poisson distribution (PD) and Marginal Poisson distribution (MPD) versus photon counts as 
given by: eqs. II17I I. ISl l and JSjl for the exponential prior pdf, panels (a)-(c); eqs. I I17II . (O and JTOj for the inverse-Gamma function 
prior pdf, panels (d)— (f). 



Equation (|18p allows one to provide conclusive posterior 
pdfs for the detected sources from combined energy bands. 
It results, as expected, that if an object is detected in mul- 
tiple bands, the resulting p(i3ij | dij)comb is larger than the 
source probability obtained analysing a single band. An ap- 
pUcation of this technique is shown in Section 18.1.11 This 
analysis can be extended to study crowded fields or blends 
by investigating source colours. 



I is the intensity of the source, i.e. the net source counts. 
(Jx,(Ty and p provide the source shape. Xs, ys is the source 
pixel position on the image. 

Position, intensity and source shape are provided max- 
imizing the likelihood function assuming constant priors for 
all parameters: 

TT ^i' 

p{xs,y^,I,cT^,cTy,p\b,d) oc I I — i^expf-Aj} 



4 SOURCE CHARACTERIZATION 

Following the estimation of the background and the identi- 
fication of sources, the sources can be parameterized. 
SPMs at different resolutions are investigated first. Sources 
are catalogued with the largest source probability reached 
in one of the SPMs. Local regions are chosen around the 
detected sources. The region size is determined by the cor- 
relation length where the maximum of the source probability 
is reached. 

Although any suitable source shape can be used, we 
model the sources by a two-dimensional Gaussian. The data 
belonging to a source detection area 'k' are given by: 



Aj=&ij+Gij V{ij}G{fc}. 



(19) 



Dij are the modelled photon counts in a pixel {ij} spoiled 
with the background counts Gij is the function which 
describes the photon counts distribution of the detected 
source: 

/ . 1 



Gij — 



27rcrxay — p 



■ exp{- 



2(1 -p2) 
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2p(- 



-)]}• 



V{ij} G {fc}. 



(20) 



where are the photon counts detected on the image. 

According to p9|) and (|20p the source fitting is executed 
on the sources for given background. This is reasonable since 
the uncertainty of the estimated background is small. No 
explicit background subtraction from the photon counts is 
needed for estimating the source parameters. 

Source position and extensions are converted from 
detector space to sky space. Source fluxes are provided 
straightforwardly. 

The rms uncertainties of the source parameters 
are estimated from the Hessian matrix, where Hij := 

— ^ Q^^^l^^^^^^ is element {ij} of the Hessian matrix and ^ 
indicates the source parameters. The square root of the diag- 
onal elements of the inverse of the Hessian matrix provides 
the la errors on the estimated parameters. 

The output catalogue includes source positions, source 
counts, source count rates, background counts, source shapes 
and their errors. The source probability and the source de- 
tection's resolution are incorporated. 

Source characterization can be extended with model se- 
lection. With the Bayesian model selection technique, the 
most suitable model describing the photon count profile 
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of the detected sources can be found. The models to be 
employed are, for instance, Gaus sian profile, King profile 
(jCavaliere fc Fusco-Femiand Il978l ). de Vaucouleurs model, 
Hubble model. Such an extention to the actual method 
would allow an improvement in the estimation of the shape 
parameters of faint and extended sources. 



5 RELIABILITY OF THE DETECTIONS 

There are several methods for reducing the number of false 
positi ves, one is to use p-values from statistical significance 
tests ('Linne mannll2003l ). P-values are used to a great extent 
in many fields of science. Unfortunately, th ey are commonly 
misin terpreted. Researches on p-values (e.g. fBerger fc Sellkd 
Il98it l showed that p-values can be a highly misleading mea- 
sure of evidence. 

In Section [5. II we provide the definition of p-values. In 
Section r5.2l we express our general view on how this problem 
can be tackled with BPT. The commonly used measure of 
statistical significance with Poisson p-values is introduced 
in Section [5.31 In Section [5.41 simulations are used for com- 
paring the Poisson p-values with the Bayesian probabilities. 

5.1 P— values 

In hypothesis testing one is interested in making inferences 
about the truth of some hypothesis -ffo, given a set of ran- 
dom variables X: X ~ f (x), where /(x) is a continuous den- 
sity and X is the actual observed values. A statistic r(X) is 
chosen to investigate compatibility of the model with the ob- 
served data X, with large values of T in dicating less compat- 
ibility (|Sellke. Bavarri fc Bergeil 120011 ). The p-value is then 
defined as: 

p = Pr(r(x)>r(x)|7?o). 

The significance level of a test is the maximum allowed 
probability, assuming Hq, that the statistic would be ob- 
served. The p-value is compared to the significance level. If 
the p-value is smaller than or equal to the significance level 
then Ho is rejected. The significance level is an arbitrary 
number between and 1, depending on the scientific field 
one is working in. However, classically a significance level of 
0.05 is accepted. Unfortunately, the significance level of 0.05 
does not indicate a particular strong evidence against Hq, 
since it just claims an even chance. 

An extensive literature dealing wi th misinterpre- 
tations about p-values exists , see e .g. 
19871) ■ 'Bcrgcr fc Dclamp advl ll 19871 ). 
19881 ), Ibelampadv fc Borgej (ll990l )rLore do (1990 1992) 



Berger fc Sellke 



Berger fc Berry 



iKass fc Raftervl l|l9S)5h . Sellke et al. (2001). iGregor^ (|2005l ) 



and references therein. 



5.2 The Bayesian viewpoint 

Since the state of knowledge is always incomplete, a hy- 
pothesis can never be proved false or true. One can only 
compute the probability of two or more competing hypothe- 
ses (or mode l s) on the basis of the only data available, see 
e.g.£regor3 dSoi). 

The Bayesian approach to hypothesis testing is concep- 
tually straightforward. Prior probabilities are assigned to 



all unknown hypotheses. Probability theory is then used to 
compute the poste rior probabil ities of the hypotheses given 
the observed data teergeJ IiggTl ). This is in contrast to stan- 
dard significance testing which does not provide such in- 
terpretation. In fact, in the classic approach the truth of a 
hypothesis can be only inferred indirectly. 

Finally, it is important to underline that the observed 
data and parameters describing the hypotheses are sub- 
ject to uncertainties which have to be estimated and en- 
coded in probability distributions. With BPT there is no 
need to distinguish between statistical (or random) and 
systematic uncertainties. Both kinds of uncertainties are 
treated as lack of knowledge. For more on the subject see 
iFischer. Dinklage fc PaschI l|2003l ). 



5.3 Significance testing with p— values 

Several measures of statistical significance with p-values 
have been developed in astrophysics. A critical compari- 
son and discussion about methods for measurin g the signifi- 
cance of a particular observation can be found in lLinnemannI 
(2003). Following Linnemann, our attention is focused on the 
Poisson probability p-value: 



PP = P(> d\b) 



by 



J'- 



(21) 



pp is the probability of finding d or more (random) events 
under a Poisson distribution with an expected number of 
events given by 6. Linnemann remarks that Poisson proba- 
bility p-value estimates lead to overestimates of the signifi- 
cance since the uncertainties on the background are ignored. 

5.4 Comparing threshold settings for source 
reliability 

In order to restrain the rates of false source detections per 
field, a threshold on probabilities is commonl y set according 
to the goal of a specific research. For instance. iFreeman et al.l 
( 2002. ) have chosen an upper limit of 1 spurious detection 
per analysed CHANDRA field. This value corresponds to 
a false-positive probability threshold of 10~^. The method 
developed by Freeman et al. is a wavelet-based source de- 
tection algorithm (wavdetect ). For the XM M-Newton 
serendipitous source catalogue, IWatson et al.l (|2003 ) have 
chosen a detection likelihood L = 10, corresponding to ~ 4o-. 
L is given by — ln(l — P), where P is the probability of 
source detection obtained by a maxim um likelihood method 
l|Cruddace. Hasinger fc SchmittI Il988l ) . The selected likeli- 
hood threshold corresponds to the detection of < 1 spurious 
source per field. In any systematic investigation to source 
detection, the threshold level is a tradeoff between the de- 
tection power and false detections rate. 

Following this idea, the Poisson probability p-value (|21|l 
is compared to the Bayesian source probability ((17} . Fig. [6] 
compares the two statistics. Panel (a) shows the relation 
between pp and (1 — Psourco) for varying background ampli- 
tudes and source intensities. Panel (6) displays the same data 
but with fixed background value and varying source inten- 
sities. These results are obtained employing the exponential 
prior. A more detailed study, incl uding the inv e rse-G amma 
function prior, can be found in iGugUelmettH l|2009l ). The 
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Figure 6. Classic iiypothesis testing versus Bayesian approach!. 



Table 2. Bayesian source probability (Psourcc) vs Poisson p-value 
(PP). 



b 

(cnts) 


A 

(cnts) 


d 

(cnts) 




1 - PP 


1 Psrc 


PP 


0.1 


1.0 


2 


0.97 


0.995 


0.03 


0.005 


1.0 


10.0 


2 


0.29 


0.74 


0.71 


0.26 


10.0 


100.0 


2 


0.01 


0.0005 


0.99 


0.9995 



This example shows the variation of Psourco, here indicated with 

Pare, and PP for detecting 2 photon counts in a pixel at three 
background levels. Psourcc represents the source probability 
of detecting 2 photon counts (d) in a pixel according to the 
Bayesian technique. (1-pp) is the cdf. It provides the probability 
of detecting 1 photon count or less in a pixel according to 
Poisson statistics. (1-Psource) is the background probability 
estimated with the Bayesian method, pp provides the probability 
of detecting 2 photon counts or more in a pixel according to 
Poisson statistics. 

squares on each plot indicate the photon counts d chosen 
for calculating (I21|l . On the abscissa, the background proba- 
bility is calculated as the complementary source probability 
provided by the Bayesian method. The value close to unity 
corresponds to a source probability, Psourcc, which goes to 
zero, instead a value of 0.1 corresponds to 90 per cent source 
probability and 0.01 to 99 per cent source probabihty. 
Each plot shows a general tendency. For a given count num- 
ber d, Psourcc and (1-pp) increase with decreasing back- 
ground intensity. However, Psourcc is more conservative. This 
is due to the dependency of Psourcc not only on the mean 
background intensity but also on the source intensity distri- 
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Figure 7. Dependency of the Bayesian source probability on the 
source intensities distributed on a field. 



bution. This dependency is expressed by the likelihood for 
the mixture model ([5)l, that plays a central role for the es- 
timation of the source probability (|17|l . An example of the 
different interpretations of source detection employing the 
two statistics is provided in Table [2] 

In Fig. [T] Psourcc for a given number of photon counts 
in a pixel versus X/b, the ratio between the mean source in- 
tensity and the background intensity, is drawn for a fixed 
background value. The abscissa is drawn on a logarithmic 
scale. For a given number of photon counts, the value of 
Psourcc varies with the source intensities expected in the as- 
tronomical image and with the background amplitude. For 2 
photon counts, Psourcc reachcs a maximum where the mean 
source intensity in the field has values in the range (1 — 5) 
counts. In this part of the curve, 2 photon counts in a pixel 
are discriminated best from the background. Away from this 
range, the source probability decreases. For small X/b values, 
Psourcc approaches 0.5 because source and background can 
not be distinguished. For large X/b values, Psourcc decreases 
since more sources with large intensities are expected rela- 
tive to small intensities. Therefore, a signal with 2 counts is 
assigned to be background photons only. 
(1 — pp) is calculated for the same values of background and 
photon counts as for Psourcc- (1 — pp) is constant, since it 
does not depend on the source intensities expected in the 
field. Its value for 2 photon counts, as seen in Table (2] is 
larger than the maximum of Psourcc- In general, (1 — pp) is 
larger than Psourcc for values of photon counts larger than 
the mean background intensity. If the values of the photon 
counts are lower than or equal to the mean background in- 
tensity, (1 — pp) is lower than the maximum of Psourcc- 

The comparison shows that it is not possible to cali- 
brate pp with Psourcc because of the intrinsic difference in 
the nature of the two statistics. In fact, Poisson p-values do 
not provide an interpretation about background or sources 
and do not include uncertainties on the background. The 
Bayesian method, instead, gives information about back- 
ground and sources and their uncertainties. 

The comparison between the two statistics reveals that 
slightly different answers are arising for the two priors of 
the source signal. When the exponential prior is employed, 
fields with large intensities are less penalized by false posi- 
tives caused by random Poisson noise than fields with source 
signal very close to the background amplitude. When the 
inverse-Gamma function prior is used, false positives detec- 
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tions depend on the cutoff parameter a. This is because the 
cutoff parameter has an effect on faint sources. The same 
behaviour is expected on false positives in source detec- 
tion. The exponential prior, instead, does not exclude small 
source intensities. Note that the choice of the source sig- 
nal prior pdf is crucial for source detection. For a reliable 
analysis the source signal prior chosen has to be as close as 
possible to the true one. 



6 SIMULATED DATA 

Artificial data are used for performance assessment of the 
BSS technique. Three simulations are analysed utilizing the 
exponential and the inverse-Gamma function prior pdfs of 
the source signal. The datasets, described in Section [6. II are 
meant to test the capabilities of the BSS method at varying 
background values. The idea is to cover different cases one 
encounters while surveying different sky regions or employ- 
ing instruments of new and old generations. In Section [6.21 
we review the outcome of our analysis on the three simu- 
lated datasets. Comments are given for each feature of the 
developed technique. In Section [6.31 we provide a summary 
on the outcome of our analysis. 



6.1 Simulations setup 

Three sets of simulated fields composed of 100 sources mod- 
elled on a constant background with added Poisson noise are 
generated. Groups of ten sources are characterized by the 
same number of photon counts but with different sizes. A 
logarithmic increment in photon counts per group is chosen 
ranging from 1 to 512. The shape of each source is charac- 
terized by a two-dimensional circular Gaussian. The source 
extensions, given by the Gaussian standard deviation, in- 
crease from 0.5 to 5.0 pixels in steps of 0.5. Sources are 
located equidistantly on a grid of 500 x 500 pixels. Initially, 
the simulated sources are located on the grid such that the 
source intensities increase on the abscissa, while the source 
extensions increase on the ordinate. Subsequently, the 100 
sources are randomly permuted on the field. A background 
is added on each simulated field with values of 0.1, 1 and 10 
counts respectively. We assume a constant exposure. 

In Fig. [8] we show one out of three datasets: the simu- 
lated data with small background. Image (a) represents the 
simulated data with added Poisson noise. The image indi- 
cated with (&) is the simulated data without Poisson noise. 
It is placed for comparison. The simulated datasets (a) and 
(6) are scaled in the range (0 — 3) photon counts-pixeP^ 
in order to enhance faint sources. The original scale of 
the simulated data with Poisson noise is (0 — 317) photon 
counts-pixel"^. The simulated data without Poisson noise 
show counts-pixel"^ in the range (0.1 — 326.0). 
The images representing the other two datasets for 6=1 and 
10 counts are similar to the one shown in Fig. (8] In these 
datasets, the number of sources to be separated from the 
background decreases with increasing background intensity. 

The cutoff parameter a is chosen to be 0.14 counts in 
the three simulated datasets. This is to show the effect of a 
when the background is smaller or larger than a. 



6.2 Results 

6.2.1 Background estimation 

For the background modelling, only four pivots located at 
the field's corners are used. This choice is driven by the pres- 
ence of a constant background. An optimization routine is 
used for maximizing the likelihood for the mixture model 
((5]). The solution of the optimization routine is the pivots 
amplitude's estimates from which the background is calcu- 
lated. 

The three setups are designed such that half of the 100 
simulated sources are characterized by ^ 16 photon counts. 
Some of these simulated sources are too faint for being de- 
tected. These sources may contribute to the background 
model. 

In Fig. [8l the estimated background maps are displayed 
when employing the exponential prior pdf (image (d)) and 
the inverse-Gamma function prior pdf (image (/)) for the 
simulated data with small background. 

The two images show that the background intensity de- 
creases slightly toward the upper left and lower right cor- 
ners of about 5 per cent. The same trend is seen also in the 
estimated backgrounds with intermediate and large values. 
Evidently, this effect is not introduced by the prior over the 
signal. Also, it is not introduced by the selected pivots po- 
sitions. If that were the case, then the same magnitude is 
expected at each image corner. Instead, this is an overall ef- 
fect induced by the simulated sources. All simulated sources 
are randomly permuted. In the upper left and lower right 
corners are located numerous faint sources. In the lower left 
corner many bright sources are clustered. The increment in 
the background intensity is due to the statistical distribu- 
tion of the sources. This explains why the same trend in 
background intensities is seen in all background models. 

When employing the exponential prior pdf, the esti- 
mated background intensities are in agreement with the sim- 
ulated background amplitudes. In the case of the inverse- 
Gamma function prior pdf, the estimated backgrounds are 
sensitive to the cutoff parameter a. When a is set larger 
than the mean background (i.e. simulated data with small 
background), the background is overestimated. The overes- 
timated background is due to the presence of source signal 
below the cutoff parameter. Hence, no source intensities be- 
low 0.14 counts are allowed. It results that the estimated 
background is 40 per cent larger than the simulated one. 
For simulated data with intermediate background, the cut- 
off parameter a is fixed to a value lower than the simulated 
background. The background is underestimated by only ~ 1 
per cent with respect to the simulated one. For simulated 
data with large background, the cutoff parameter a is much 
lower than the simulated mean background value. The esti- 
mated background is in agreement with the simulated one. 

The background uncertainties are quite small compared 
with the background itself, on the order of few a percent. 
This effect holds because the background is estimated on the 
full field. However, the errors increase where the estimates 
deviate from the simulated background. In addition, when 
applying the inverse-Gamma function prior pdf, the errors 
are larger than those found utilizing the exponential prior 
pdf. 
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Figure 8. Panels (a) and (b), simulated data with small background: image with Poisson noise, image without Poisson noise, respectively. 
Panels (c) and (d), results with exponential prior pdi: SPM with 3 pixels resolution, background map, respectively. Panels (e) and (/), 
results with inverse-Gamma function prior pdf. 



6.2.2 Hyper-parameter estimation 

In Fig. [9] the contour plot in (A, /3) parameter space for 
the joint probabihty distribution is shown for the hyper- 
parameters evaluated from the simulated data with small 
background. The contour levels indicate the credible re- 
gions. The values of the estimated hyper-parameters are: 
15 = (99.2 ± 0.03) per cent, A = (3.68 ± 0.1) counts. The 
estimated /3 value provides the information that only 0.8 
per cent of the pixels in the field contains sources. A similar 
answer is found with the other simulated data: the (5 value 
increases slightly at increasing background amplitudes. A, 
instead, provides the mean source intensity in the field. The 



estimated value of A increases with increasing background 
amplitudes because small intensities are assigned to be back- 
ground. 

When employing the inverse-Gamma function prior 
pdf, the hyper-parameter a is found with the smaller value 
in the simulated data with small background. The largest 
value of a is found in the simulated data with intermedi- 
ate background. Large values of a indicates that more faint 
sources and less bright sources are expected in the field 
(Fig. [TJ. These results do not contradict our expectations 
on the hyper-parameter estimates, since the cutoff parame- 
ter selects the source signal distribution at the faint end. 
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Figure 9. Contour plot of the posterior pdf for the hyper- 
parameters, p{\, I3\D), estimated from the simulated field with 
small background using the exponential prior pdf. 



6.2.3 The components of the mixture model 

In Fig. 1101 the Poisson and the marginal Poisson distribu- 
tions multiplied with their prior knowledge on the model 
are plotted over the normalized histogram of each simulated 
dataset. The likelihoods are drawn for the exponential and 
inverse-Gamma function prior pdfs. The values of the esti- 
mated hyper-parameters are shown in each image. 

The exponential prior pdf is plotted over the histogram 
with a continuous line. The inverse-Gamma function prior 
pdf, instead, is plotted with a dashed line. The simulated 
data are neither distributed exponentially nor as an inverse- 
Gamma function. Hence, the prior pdfs of the source signal 
are not expected to fit the data exactly. 

The marginal Poisson distribution weighted with (1 — /?) 
is drawn with a dash-dot line when employing the exponen- 
tial prior and with long dashes line in the case of the inverse- 
Gamma function prior. The Poisson distribution (dotted 
line) is weighted with /3 for the exponential and the inverse- 
Gamma function prior pdfs. 

The intersection between the Poisson pdf and the 
marginal Poisson pdf indicates the source detection sensi- 
tivity. When employing the exponential prior in the sim- 
ulated data with small background (panel (a)), the expo- 
nential prior enables the detection of fainter sources than 
the inverse-Gamma function prior. This is expected since 
the cutoff parameter occurs at a value larger than the sim- 
ulated mean background. Considering the simulated data 
with intermediate background (panel (6)), the detection is 
more sensitive to faint sources when employing the inverse- 
Gamma function prior compared to the exponential prior. 
In fact, the cutoff parameter allows to describe as source 
signal part of the simulated background amplitude. For the 
simulated data with large background (panel (c)), the same 
sensitivity in source detection is expected when employing 
the two priors over the signal distribution. 



6.2.^ Source probability maps 

The box filter method with cell shape of a circle is used in 
the three simulations for the multi-resolution analysis. Ex- 
amples of SPMs are shown in Fig. [S] for the simulated data 
with small background. Images (c) and (e) are obtained em- 
ploying the exponential and the inverse-Gamma function 
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Figure 10. Normalized histograms of the simulated data with 
small (panel (a)), intermediate (panel (6)) and large (panel (c)) 
backgrounds are displayed. The exponential prior (Epr), inverse- 
Gamma function prior (IGpr) and related marginal Poisson 
(MPD) and Poisson (PD) distributions are plotted over the data. 
The ordinates are in logarithmic scale. 



prior pdfs, respectively. These images represent the prob- 
ability of having source contributions in pixel cells with a 
resolution of 1.5 pixels. At this resolution a pixel cell is com- 
posed by 9 pixels. A pixel cell with a correlation radius of 
1.5 pixels is drawn in the lower right corner of image (c) 
(Fig. [SI). It is indicated with an arrow. 

The multi-resolution technique provides an analysis of 
source probabilities variation and of source features of the 
detected sources. In Fig. 1111 we display the photon count 
image and the SPMs zoomed in a bright extended source. 
This source is detected with the largest source probability 
(~ 1) at 2.5 pixels resolution. At this resolution the source 
is detected as one unique object, as given by the simulation. 
At larger resolutions, instead, the source is dissociated in 
small parts as seen in the photon count image. This indi- 
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Figure 11. The upper left image is a zoom in of the photon count image (panel (a), Fig. [S]! on a simulated source located at 
(x,y)=(360,270). The width of the photon count image is 65 pixels. The following images are SPMs at decreasing resolutions. The 
correlation length of each SPM is written on the lower right corner of each image in pixel units. 



cates that the source counts are not distributed uniformly. 
In an astronomical observation, more information is required 
in order to understand the nature of such sources. Secondly, 
the maximum in source probability is reached at a correla- 
tion length that is smaller than the source size. This is due 
to the source brightness relatively to the small background 
value. Within the range of resolutions studied, the source 
probability is constant at correlation lengths larger than 2.5 
pixels. 

This example shows that the multi-resolution technique 
combined with the BSS method is particularly appropriate 
for the search of extended and non-symmetrical sources. 

6.2.5 Comparison between estimated and simulated source 
parameters 

Source parameters and their uncertainties are derived as de- 
scribed in Section U Sources are catalogued when a prob- 
ability larger than 50 per cent is reached at least in one of 
the SPMs. This threshold is chosen within these simulated 
datasets in order to clarify the interpretations provided in 
Tabled 

The parameters of bright sources are precisely esti- 
mated. In Fig. [TT] we give an example of detection of 
a bright extended source employing SPMs in the multi- 
resolution technique. The true parameters of this source 
are: 128 photon counts, (ax,o-y)=(5,5), (s,y) = (360,270). The 
estimated parameters of this source are: (129.79 ± 23.70) 
net source counts, (o-x,(Ty) = (4.92 ± 0.67,4.96 ± 0.71) pixels, 
(a;,y) = (359.57±0.92,269.29±0.98) pixels. Instead, the effect 
of background fluctuations on faint source estimates can be 
quite pronounced. 

In Fig. [T21 we provide an example of the estimated 
source positions and extensions on the simulated data with 
small background, utilizing the exponential prior pdf. The 
errors on the estimated parameters are not considered in this 
plot. Some of the detected faint sources look uncentered and 
distorted. Four false positives in source detection are found. 



Table 3. Source detections on simulated data employing different 
prior pdfs of the source signal. 
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prior pdf 
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The background value of the simulated data are reported in 
counts-pixel"^. Epr and IGpr have same meaning as given in 
Fig. 1101 The term true detect provides the number of detected 
sources matched with the simulations. They are catalogued if 
their probabilities are larger than 50 per cent. The number of 
false positives in source detection are listed considering 50, 90 
and 99 per cent probability thresholds. 

The simulated data without Poisson noise with the simu- 
lated source shapes superimposed are shown for comparison. 

Table |3] reports the number of detected sources for each 
simulation. Different columns are used for accounting true 
detections and false positives separately. The number of de- 
tected sources employing the inverse-Gamma function prior 
is larger with respect to the exponential prior case only when 
the cutoff parameter is set lower than the mean background 
amplitude. 

In Fig. 1131 the estimated source counts are related to 
their resolution of source detection (panels (a) and (c)) and 
their source probability (panels (6) and (d)) for the sim- 
ulated data with small background. These are log-linear 
plots. The results utilizing the exponential and the inverse- 
Gamma function prior pdfs are shown in panels (a) and (6) 
and in panels (c) and (d), respectively. The asterisks indicate 
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Figure 12. Panels (a) — (b): simulated data with small background with and without Poisson noise, respectively, and scaled in the range 
(0 — 3) counts. The catalogue obtained employing the exponential prior pdf (a) and the shapes of the region of the simulated sources (b) 
are superposed. 
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Figure 13. Results on simulated data with small background employing the exponential prior pdf (panels (a) and (6)) and the inverse- 
Gamma function prior pdf (panels (c) and {d)). Panels (a) and (c): correlation length in pixel units versus the net source counts. Panels 
(fe) and (d): source probability versus net source counts. Sources matched with the simulated input catalogue are indicated with an 
asterisk. A square indicates false positives in source detection. 



true sources, while the squares show spurious detections. 
The effect of the cutoff parameter a on source detection is 
visible in panel (c). In this example, the value of a is chosen 
larger than the simulated background amplitude. Hence, the 
inverse-Gamma function prior pdf does not allow to detect 
sources as faint as the exponential prior. 
The plots on panels (a) and (c) denote the source detec- 
tion technique at several correlation lengths. In both plots, 
a line is drawn only with the purpose to guide the eye. Left 



to the line, sources can not be detected. Right to the line, 
sources can be found. Very faint sources with few photon 
counts 10) are resolved by small correlation lengths. In 
the range (10 — 100) net source counts, sources are detected 
at decreasing resolutions. Bright sources do not require large 
correlation lengths for being detected. 

The plots in panels (fe) and (d) of Fig. llSl highlight that most 
of the detections occurs at probabilities larger than 99 per 
cent. At this probability value and larger, sources matched 
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Figure 14. Relation between the simulated sources with small background and the measured sources as in output from the processing 
with the developed method employing the exponential prior pdf. Plots (a) and (b) show the comparison of the measured source positions 
with the simulated input positions on the X axis and on the Y axis. Plot (c) displays the relation of the measured source photon counts 
versus the simulated intensities. A comparison of the estimated source extensions (Tx and ay with the input source extensions are displayed 
in the plots (d) and (e). The errors estimated for the source parameters are overplotted. The error bars on the estimated values denote 
the 68 per cent confidence limit of the corresponding posterior distribution. The lower panel in each image shows the fractional residuals. 



with the simulated one are strongly separated from false 
positives. 

In Fig. [14] we show the relation between the simulated 
and the estimated source parameters. Good estimates in 
source parameters are achieved. The estimated net source 
counts errors and extension errors can be large for faint 
sources. 

In Fig. [12] a summary on the analysis of all the detected 
sources employing the exponential prior on the three sim- 
ulated datasets is presented. The plot in panel (a) shows 
the difference between estimated and simulated net source 
counts, normalized with the estimated errors, versus the 
source probability of the merged data. This plot does not 
contain the information provided by false positives in source 
detection. 87 per cent of all detections occurred with prob- 
ability larger than 99 per cent. The image in panel (b) is 
a semi-log plot of the normalized difference between mea- 
sured and simulated net source counts versus the simulated 
net source counts of 98.5 per cent true sources detected in 
the three simulations. The values of two sources, detected in 
the simulated data with large background, are outside the 
selected y range. These two detections are included in the 
analysis of verification with existing algorithms (Section [7]). 
The residuals are normally distributed, as expected. They 
are located symmetrically around zero. At the faint end, 
the results are only limited by the small number of simu- 
lated faint detectable sources. Faint and bright sources are 
equally well detected. 



6.2.6 False positives 

Until now we have discussed the detections that have coun- 
terparts with the simulated data. We consider now the de- 
tection of false positives. 

In Table [S] the number of detected false positives are 
listed for each simulation. At 50 per cent probability thresh- 
old more false positives are found with the inverse-Gamma 
function prior compared to the exponential one. At a 90 
per cent source probability threshold, the analyses with the 
two prior pdfs provide similar results. True detections are 
strongly separated from statistical fluctuations for source 
probability values larger than 99 per cent. 

When employing the inverse-Gamma function prior 
pdf, the number of false positives is sensitive to the cut- 
off parameter. Less false positives are found when a is set 
larger than the background, because it reduces the number 
of detectable faint sources. It may be worth noting that even 
when the cutoff parameter is set larger than the background, 
a probability threshold ^ 90 per cent has to be considered 
(see Fig. [13] for more details) . 

False positives in source detections show large errors 
in their estimated parameters. The source probability varia- 
tion and the source features analyses in the multi-resolution 
technique provide hints of ambiguous detections. However, 
as all methods, the BSS approach is limited by statistics. 
Spurious detections can not be ruled out. 
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Figure 15. Merged information from three simulated datascts of 
the detected sources employing the exponential prior pdf. Panel 
(a) , difference between estimated and simulated net source counts 
normalized by the errors on the estimated net source counts ver- 
sus source probability. Panel (6), difference between estimated 
and simulated net source counts normalized by the errors on the 
estimated net source counts versus the simulated counts. 



6.2. 7 Choice of the prior pdf of the source signal 

The big difference between the two prior pdfs of the source 
signal follows on from one prior pdf having one parameter 
and the other pdf having two. 

The parameter A, indicating the mean intensity in an 
astronomical image, introduced with the exponential prior 
pdf is estimated from the data. 

The parameter a, that is the shape parameter of the 
power-law, given by the inverse-Gamma function prior pdf 
is estimated from the data. Instead the cutoff parameter 
a is selected to a small value such that the inverse-Gamma 
function prior pdf behaves as a power-law. Astronomical im- 
ages can be characterized by a small background. It results 
that a can be chosen from a number of alternatives, rang- 
ing from values that are above or below the background 
amplitude. The choice of a implies a selection on the de- 
tectable sources: sources whose intensity is lower than a are 
not detected; sources close to the background amplitude are 
detected when a is set below the background amplitude. 

On real data much more prior information for the cut- 
off parameter is needed. The inverse-Gamma function prior 
pdf can be employed if a mean value of the background am- 
plitude is already known from previous analyses. 

The exponential prior pdf is preferable over the inverse- 
Gamma function prior pdf, since no predefined values are 
incorporated. This is also supported by the results obtained 



with the simulated data. However, the inverse-Gamma func- 
tion prior is a more suited model to fit the data and it has 
potentials for improving the detections of faint objects. 

One way to improve the knowledge acquired with the 
inverse-Gamma function prior is by the estimation of the 
cutoff parameter from the data. This change in our algo- 
rithm is not straightforward and MCMC algorithms have to 
be employed. 

6.3 Sumrriciry 

Simulated data are employed to assess the properties of the 
BSS technique. The estimated background and source prob- 
abilities depend on the prior information chosen. A success- 
ful separation of background and sources must depend on 
the criteria which define background. Structures beyond the 
defined properties of the background model are, therefore, 
assigned to be sources. There is no sensible background- 
source separation without defining a model for the object 
background. Additionally, prior information on source inten- 
sity distributions helps to sort data, which are marginally 
consistent with the background model, into background or 
source. Therefore, prior knowledge on the background model 
as well as on the source intensity distribution function is cru- 
cial for successful background-source separation. 

For the background model a two-dimensional TPS rep- 
resentation was chosen. It is flexible enough to reconstruct 
any spatial structure in the background rate distribution. 
The parameters are the number, position and amplitudes of 
the spline supporting points. Any other background model 
capable to quantify structures which should be assigned to 
background can be used as well. 

For the prior distribution of the source intensities the 
exponential and the inverse-Gamma function are used as 
illustrations. For both distributions the source probability 
can be given analytically. The hyper parameters of both 
distributions can either be chosen in advance to describe 
known source intensity properties or can be estimated from 
the data. If they are estimated from the data simultaneously 
with the background parameters, properties of the source 
intensity distribution can be derived, but at the expense of 
larger estimation uncertainties. It is important to note that 
the performance of the BSS method increases with the qual- 
ity of prior information employed for the source intensity 
distribution. The prior distribution of the source intensities 
determines the general behaviour of the sources in the field 
of view and the hyper-parameters are useful for fine-tuning. 

The aim of detecting faint sources competes with the 
omnipresent detection of false-positives. The suppression of 
false-positives depends both on the expedient choice of prior 
information and on the level of detection probability ac- 
cepted for source identification. Compared to, e.g., p- values 
the BSS technique is rather conservative in estimating source 
probabilities. Therefore, a probability threshold of 99 per 
cent is mostly effective to suppress false positives. 

The estimated background rates are consistent with the 
simulated ones. Crowded areas with regions of marginally 
detectable sources might increase the background rate ac- 
cordingly. 

The SPMs at different correlation lengths are an impor- 
tant feature of the technique. The multi-resolution analysis 
allows one to detect fine structures of the sources. 
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Table 4. Number of detected sources employing WAVDETECT on 
three simulated datasets. 



simulated 


true detect 


false 


data 




positives 


0.1 


56 


4 


1.0 


37 


1 


10.0 


23 


1 



The explanation of these columns can be found in Table |3] 

The source parameters are well determined. Their resid- 
uals are normally distributed. We will show in Section [7] 
that the BSS technique performs better than frequently used 
techniques. Naturally, the estimation uncertainties of pa- 
rameters for faint sources are large due to the propagation 
of the background uncertainty. 



7 VERIFICATION WITH EXISTING 
ALGORITHMS 

In the X-ray regime, the sliding window technique and the 
WT techniques are widely used. However, the WT has been 
shown to perform better than the sliding window technique 
for source detection. The WT improvement in source detec- 
tion with respect to the sliding window techni que is inversely 
prop ortional to the background amplitude (jPreeman et al.l 
l2002f ). The WT has also other favourable aspects for being 
compared with the BSS method developed in our work. The 
WT allows for the search of faint extended sources. The BSS 
method with the multi-resolution technique is close to the 
WT method. 

Between all the available software e mploying WT, we 
choose WAVDETECT (iFreeman et al.|[2002l ). part of the freely 
available CIAO software package. We use version 3.4. 

WAVDETECT is a powerful and flexible software package. 
It has been developed for a generic detector. It is applica- 
ble to data with small background. The algorithm includes 
satellite's exposure variations. It estimates the local back- 
ground count amplitude in each image pixel and it provides 
a background map. 

We utilize wavdetect on our simulated data described 
in Section |6] The threshold setting for the signiflcance 
('sigthresh') is chosen equal to 4.0 x 10~^, in order to de- 
tect on the average 1 spurious source per image. The 'scale' 
sizes are chosen with a logarithmic increment from 2 to 
64. Tests have been made changing the levels of these pa- 
rameters, assuring us that the selected values provide good 
performance of this WT technique. In Table |4l we report 
the number of detected sources per simulated field, sepa- 
rating the sources matched with the simulated one {true 
detect) to the false positives in source detection found em- 
ploying the above mentioned threshold setting. The three 
simulated datasets are distinguished by their background 
values (counts). The simulated background values are re- 
ported in column simulated data. These results are com- 
pared with the ones obtained with the BSS algorithm when 
employing the exponential prior pdf as shown in Table [3] 



The BSS technique finds all sources detected by wavde- 
tect. In the simulated data characterized by small back- 
ground, the two algorithms find the same number of false 
positives. The BSS algorithm reveals 8 per cent more true 
detections than wavdetect. These sources are character- 
ized by counts in the range (4 — 8). Hence, the BSS tech- 
nique performs better than wavdetect in the low number 
counts regime. The performance of the two techniques on 
the datasets with intermediate and large background values 
is similar. 

The explanation for these results is enclosed in the back- 
ground estimate. The wavdetect estimates of background 
values are similar to the results obtained with the BSS tech- 
nique in the intermediate and large background datasets. 
Though, the backgrounds provided by wavdetect present 
rings due to the Mexican Hat function employed as a filter 
on the image data. In the simulated data with small back- 
ground, the WAVDETECT background model is on the average 
larger than the one estimated with the BSS method. The 
plots in Fig. [16] support these conclusions (semi-log plots) . 
The image in panel (a) shows that wavdetect fluxes are 
underestimated for ~ 20 per cent of all detected sources. In 
addition, wavdetect sensitivity for source detection is lim- 
ited to 16 counts per source within these simulated datasets. 
In Fig. [15] panel (b) presents the sensitivity achieved by the 
BSS method. The plot in panel (6) of Fig. [16] displays a 
relation between the sources detected by wavdetect (ordi- 
nate) and by BSS (abscissa), both matched with the simu- 
lated data. Most of the wavdetect underestimated sources 
are coming from the simulated data with small background. 
BSS presents only two sources underestimated and detected 
in the simulated data with large background. By chance, 
the triangle located at (—8,-2) indicates the detection of 
two sources. Both sources where simulated with 256 source 
counts and a circular extension of 4 pixels one, 5 pixels the 
other. The estimated source positions are also improved with 
BSS (Fig. [13 semi-log plots). 

The residuals provided by the BSS technique are a factor of 
10 smaller than the ones from wavdetect. Their estimates 
have many outliers. Our estimates are normally distributed. 

Though the comparison between the two detection 
methods is not yet carried out on real data, these results 
are encouraging. The BSS method detects at least as many 
sources as wavdetect. The simulations prove that the de- 
veloped Bayesian technique ameliorates the detections in the 
low count regime. The BSS estimated positions and counts 
are improved. Finally, we expect that the BSS technique 
will reflne wavdetect sensitivity on real data, because the 
BSS technique is designed for modelling highly and slowly 
varying backgrounds taking into account instrumental struc- 
tures. 



8 APPLICATION TO OBSERVATIONAL DATA 

In this Section we employ RASS data to show the effective- 
ness of our method at large varying background, exposure 
values and where exposure nonuniformities are consistently 
taken into account. We present results obtained combining 
SPMs at different energy bands. We demonstrate that our 
algorithm can detect sources independently of their shape, 
sources at the field's edge and sources overlapping a diffuse 



22 F. Guglielmetti et al. 



D 
> 
< 



■10 - 



-15 - 



-20 t 



10 100 
simulated counts 
(a) 



D 
> 
< 


























- ob=0.1cts 
□ b=1.0cts 
: Ab=10.cts 









-10 -8 -6 -4-2 2 
residuals/BSS error counts 
(b) 



Figure 16. Panel (a): normalized difference between WAVDETECT estimated source counts versus simulated counts of all the detected 
sources in the three simulated datasets. Panel (6): comparison between the normalized differences of WAVDETECT source counts and 
simulated counts, on the ordinate, and the normalized difference of BSS source counts and simulated counts, on the abscissa. Sources 
detected from the simulated data with small background are indicated with a diamond. A square is used to highlight sources detected 
in the dataset with intermediate background. The detected sources in the simulated data with large background are indicated with a 
triangle. Dashed lines are drawn as a borderline of the ±3(T detection. The zero line is indicated with a continuous linestyle. 
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Figure 17. Same explanation as for panel (b) in Fig. 1161 but comparing the estimated source positions: x axis (panel (a)), y axis (panel 
{&))■ 



emission. Last, the BSS technique provides evidence for ce- 
lestial sources not previously catalogued by any detection 
technique in the X-ray regime. 



8.1 ROSAT PSPC in Survey Mode data 

We apply the BSS algorithm on a data sample coming 
from the Position Sensitive Proportional Counter (PSPC) 
on board of the ROSAT satellite in survey mode (0.1 — 2.4 
keV). T^OSAr satellite provided the only all-sky survey re- 
alized using an imaging X-ray telescope (Voges et al. 1999, 
2000). 

RASS data supply a unique map of the sky in the X-ray 
regime. The map of the sky was divided in 1, 378 fields each 
of 6.4° X 6.4°, corresponding to 512 x 512 pixels. In addition, 
the satellite's exposure time can vary between about 0.4 to 
40 ks and some parts of the sky are without observations 
due to the satellite's crossing of the auroral zones and of the 
South Atlantic Anomaly. Therefore, RASS data provide a 
wide range of possibilities for testing the BSS algorithm. 

We present the results obtained analysing three ROSAT 
fields, whose IDs are: RS930625n00, RS932209n00 and 
RS932518n00. 



8.1.1 RS930625n00 

This ROSAT field is located at a = 17''49™5^ 5 = 
-|-61°52'30". It is characterized by large variations of the 
satellite's exposure ranging from 1.7 to 13.5 ks. The count 
rate image of RS930625n00 in the broad energy band (0.1 — 

2.4 keV) is displayed in panel (a) of Fig. 1181 The count rates 
range over (0. — 0.11) photon count-s~^ -pixel" ^. The image is 
here scaled in the range (0. — 0.005) photon count-s~^ -pixel" ^ 
in order to enhance the sources. 

Our results employing the BSS algorithm are shown in 
Fig.[T8]with a SPM, the TPS map and the background map. 

In panel (b), the SPM is obtained combining statisti- 
cally the soft (0.1 - 0.4 keV) and the hard (0.5 - 2.4 keV) 
energy bands. The SPM here displayed is obtained employ- 
ing the exponential prior pdf and accounts for the width of 
the ROSAT instrumental PSF. The information of pixels is 
combined with the box filter method with the cell shape of 
a circle. This image corresponds to a correlation length of 

1.5 arcmin. Sources are identified in terms of probabilities. 
The image is in linear scale. 

In panel (c), the TPS map is modelled from 
RS930625n00 field, broad energy band. The TPS mod- 
els the background rate. Only 25 support points are 
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Figure 18. Analysis of RS930625n00, ROSAT PSPC in Survey Mode field. 



used. The background rate ranges (0.0005 — 0.001) photon 
count-s" -pixel" . The contours are superposed for enhanc- 
ing the features relative to the modelled background rate. 
The innermost and the outermost contours indicate a level 
of 0.0005 and 0.001 photon count-s~^ -pixel" , respectively. 

The corresponding background map estimated from the 
selected ROSAT field is displayed in panel (d). Its values are 
in the range (1.17—8.53) photon count - pixel" ^. The image is 
scaled in the range (1.17 — 6.68) photon count-pixel"'^. The 
contour levels close to black and light gray delineate 6.0 
and 2.0 photon count-pixel"^, respectively. The background 



map shows the prominent variation due to the heterogeneous 
satellite exposure time. 

The lower row of Fig. [TH] shows the background rate 
(panel (e)) and the background amplitude (panel (/)) as ob- 
tained analysing the broad energy band with the Standard 
Analysis Software System (SASS). SASS is the detection 
method utilized for the realization of RASS. SASS combines 
the sliding window technique with the maximum likelihood 
PSF fitting method for source detection and characteriza- 
tion, respectively. The background rate image is in the range 
(0.0005-0.00097) photon count-s"^-pixel"^ The innermost 
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Figure 19. POSS-II red plate with overlaid ROSAT contours 
corresponding to 3, 4, 4.5 and 5 cr above the local background. 
In the image centre is located SDSS J172459.31+642424.0, a low 
redshift QSO. 

and the outermost contour levels delineate 0.0006 and 0.0009 
photon count-s^^-pixel~^, respectively. The background am- 
plitude image estimated by SASS has values in the range 
(1.11 — 8.27) photon count-pixel"^. The image is scaled as 
the one obtained from the BSS technique. The contour levels 
indicate from 2.0 to 6.0 photon count-pixel"^. 

The BSS technique allows for more flexibility in the 
background model and the edges are more stable than the 
ones obtained with SASS. Hence, celestial objects located 
on the edges are not lost during source detection. 

In Fig. [19] we show an example of a QSO detection 
with the BSS algorithm analysing the field RS930625n00 
in the hard (0.5 — 2.4 keV) band image. No counterparts 
have been found in the ROSAT Bright and Faint source 
catalogues, Voges et al. (1999, 2000). A counterpart is 
found with the Sloan Digital Sky Surve y. This QSO is 
catalogued as SPSS J172459.31+642424.0 jAbazaiian et"al] 
|2004 iHao et al.ll2005l . IVanden Berk et al.ll2006l ). It is located 
at the image centre. The image covers a held of view of ~ 5 
arcmin at the side. The 7? 05^ T image, as shown in Fig. 1181 
has 45 arcsec-pixel"^ resolution. The QSO optical position 
as given by the SDSS is 12 arcsec far from the BSS position, 
i.e. it is within the pixel resolution of the ROSAT data. The 
estimated source count rate with the BSS algorithm for this 
source is (0.0068 ±0.0023) photon count-s"^ Its source de- 
tection probability is 0.999. This object is located close to 
the north-west corner of the field RS930625n00 (~ 45 ar- 
cmin). In this region the SASS background intensity is 7 per 
cent lower than the BSS background estimate. 



8.1.2 RS932209n00 

The detection capabilities of our Bayesian approach on im- 
ages with exposure non-uniformities is presented in Figs. 1201 
andim The analysed ROSAT is located at a = 3''3T", 
5 = -28°07'08". 

In Fig. [201 panel (a), the soft band image (0.1-0.4 keV) 
is shown. The image accounts for photon count-pixel"^ in 
the range (0—9). The image is scaled in the range (0—2) pho- 
ton count-pixel"^. In panel (6), an SPM as in output from 



the BSS method is displayed. The inverse-Gamma function 
prior pdf is used for source detection and background esti- 
mation. This SPM is obtained with a correlation length of 
270 arcsec. The box filter method with the cell shape of a 
circle is employed. The SPM shows source detections. 

In Fig. 1211 panel (a), the satellite's exposure is located. 
The exposure values are in the range (0. — 0.393) ks. In 
panel (6) we show the background map estimated with the 
developed BSS technique. Only 16 pivots equidistantly dis- 
tributed along the field are employed. The background am- 
plitude ranges from to 0.318 photon count - pixel" ^. The 
background is estimated with null value where no satel- 
lite's exposure information is provided. The estimated back- 
ground map is similar to the exposure map because the con- 
tribution from the cosmic background is very small. 

No artefacts, due to exposure non-uniformities, occur 
in the background map nor in the SPMs. 

8.1.3 RS932518n00 

We employ the ROSAT field id RS932518n00 in the hard 
energy band (0.5 — 2.4 keV) to show the capabilities 
of our technique on data with large background varia- 
tions. This field is characterized by a complex structure 
of hot gas caused primarily by the emissions of two galac- 
tic X-r&j sources: (1) Vela supernova remnant (SNR) lo- 
cated at a = 8''35'"20.6^ 5 = -45°10'35"; (2) SNR 
RX JQ852.0 - 4622 located at a = 8^52"* S = -46° 22', see 
refs. * Aschenbach. Egger fc Trumped l|l995h and I AschenbachI 
( 199^^ 

In Fig. 1221 panel (a), we display the photon count image 
of RS932518n00 field in the hard energy band. The image 
presents values in the range (0— 136) photon count-pixel"^. 
It is scaled in the range (0 — 15) photon count-pixel"^ in or- 
der to enhance the X-ray emissions. The satellite's exposure 
map is displayed in Fig. 1231 panel (a). The exposure map 
shows variations in the range (0.5 — 0.8) ks. 

The analysis of RS932518n00 in the hard energy band 
with the BSS algorithm is shown in Figs. 1221 panels (6)-(d), 
and 1231 panel (6) . For background estimation and source 
detection the exponential prior pdf is employed. 

In Fig. [221 panels {b)-{d), three SPMs are displayed. 
The correlation lengths used for their realization is written 
on each image. The information of neighbouring pixels is put 
together with the Gaussian weighting method. These SPMs 
are in linear scale. 

The background map fFig. 1221 panel (6)) is estimated 
utilizing 36 pivots equidistantly spaced along the field. The 
estimated background values range from ~ to 4.1 photon 
count-pixel"'^. The heterogeneous background is recovered. 

The BSS method combined with the multi-resolution 
analysis allows one to identify the detection of point-like 
sources on top of the diffuse emission and of extended X-ra,y 
features, such as the Pencil Nebula located at a = 9''0.2"', 
S = -45° 57' and the SNR i?XJ0852.0 - 4622. The BSS 
algorithm detects about 50 objects on this T^OS^r field. 

In Table [SI we present part of this catalogue. Columns 
RA, Dec, sctr, ctx, (Jy and p are the estimated positions, 
source count rates and source parameters as in output from 
the developed BSS technique. Each listed object is detected 
with a probability larger then 0.999. The column indicated 
with matched ID corresponds to an ID given by catalogues 
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RASS 3/22/09 (0.1-0.4 keV) BSS: SPM (270 arcsec) 
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Figure 20. Panel (a): RO SAT all-shy survey photon count image of the field id RS932209n00 in the soft energy band. Panel (6): source 
probability map as obtained with the developed BSS method. 




Figure 21. Panel (a): exposure map of the ROSAT field id RS932209n00 (Fig. I20|l. The z axis indicates the exposure time in s pixcl"^. 
Panel (b): estimated background amplitude with the BSS technique of RS932209n00 field in the soft energy band. The z axis shows the 
estimated background count-pixel"^. 



created analysing ROSAT data and whose position matches 
with the one obtained with the BSS algorithm. During this 
selection, we give priority to IDs provided by RASS cata- 
logues when available. Catalogues IDs highlighted with an 
asterisk are coming from a point source catalogue generated 
from all ROSAT PSPC pointed observations. 



9 CONCLUSIONS 

We have presented a new statistical method for source detec- 
tion, background estimation and source characterization in 
the Poisson regime. With this technique we have elaborated 
a very general and powerful Bayesian method applicable to 
images coming from any count detector. It is particularly 
suitable for the search of faint and extended sources in im- 
ages coming from instruments of a new generation. We apply 
our technique to X-ray data. 



The developed technique does not loose information 
from the data and preserves the statistics. The BSS algo- 
rithm provides a comprehensive error analysis where uncer- 
tainties are properly included in the physical model. 

The background and its errors are estimated on the 
complete field with a two-dimensional spline and the knowl- 
edge of the experimental data and the exposure map. The 
background model is well-defined along the complete field, 
so that field edges and instrumental effects are correctly 
handled and sources located at the edges or gaps are not 
penalized. 

The BSS method does not assume anything about the 
source size or shape for an object to be detected. Sources are 
separated from the background employing BPT combined 
with the mixture model technique. Point-like and extended 
sources are detected independent from their morphology and 
the kind of background. Consequently, the BSS technique 
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Figure 22. Panel (a): ROSAT all-sky survey photon count image of the field id RS932518n00 in the hard energy band. The analysis of 
this ROSAT iie\d with the developed BSS method is shown in panels (6) — (d) with SPMs obtained with increasing correlation lengths. 




Figure 23. Panel (a): exposure map of the ROSAT field id RS932518n00 (Fig. I22I I. The z axis indicates the exposure time in 
seconds-pixel"^. Panel (6): estimated background amplitude with the BSS technique of RS932518n00 field in the hard energy band. 
The z axis shows the background counts-pixel"^. 
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Table 5. BSS sample catalogue from the analysis of RS932518n00 field in the hard energy band (0.5 — 2.4 keV) 



RA 


Dec 


sctr 




ay 


P 


matched ID 






(count-s"-'- ) 


(pixel) 


(pixel) 







(1) 08''35'"20.4'' -45°10'36.37" 2.7259±0.0745 1.91±0.05 1.93±0.05 -0.05±0.04 IRXS J083520.6-451035(t) 

(2) 08''37'"24.62» -42°53'55.75" 0.6236±0.0486 1.49±0.17 1.56±0.12 -0.32±0.10 IWGA J0837.3-4253(*) 

(3) 08''35'"55.85'' -43°11'11.29" 0.7099±0.0524 1.55±0.12 1.73±0.15 0.06±0.11 IWGA J0835.9-431l(*) 

(4) 08'"57"'44.81^ -41°51'51.52" 0.2478±0.0253 3.20±0.33 3.60±0.38 0.07±0.14 IRXS J085748.6-41510l(t) 



References; (f) IVogea et all lll999l); (*) IWhite. Giommi fc Angelinil lll994l) 



can be applied to large data volumes, e.g. surveys from X- 
ray missions. 

We developed the BSS technique utilizing two prior pdfs 
of the source signal. We proved that the prior pdfs of the 
source signal allow one to select what has to be described as 
source or as background. 

The SPMs are robust features of our technique. They 
allow for the analysis of faint and extended objects. SPMs 
obtained at different energy bands can be combined proba- 
bilistically providing conclusive likelihoods for each detected 
source. 

The BSS photometry is robust. The errors on the esti- 
mated source parameters are normally distributed. 

We demonstrated that our technique is capable of 
coping with spatial exposure non-uniformities, large back- 
ground variations and of detecting sources embedded on a 
diffuse emission. We expect to handle consistently the het- 
erogeneities present in astronomical images with CCD pat- 
terns, images superposed and mosaic of images. 

Finally, the verification procedure with existing algo- 
rithms through simulations demonstrates that the BSS tech- 
nique improves the detection of sources with extended low 
surface brightness. This is supported by the applications on 
real data. The BSS method has good potentials for improv- 
ing, first, the count rate of detected sources and, second, the 
sensitivity reached by other techniques also on real data. 
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10 ACRONYMS 

The following is a list of acronyms frequently appearing in 
the text. 



BPT 


Bayesian Probability Theory 


BSS 


Background-Source Separation 


CCD 


Charge-Coupled Device 


cdf 


cumulative distribution function 


FWHM 


Full- Width at Half Maximum 


MCMC 


Markov-Chain Monte Carlo 


pdf 


probability density function 


PSF 


Point-Spread Function 


PSPC 


Position Sensitive Proportional Counter 


RASS 


ROSAT All-Sky Survey 


SASS 


Standard Analysis Software System 


SPM 


Source Probability Map 


SNR 


supernova remnant 


TPS 


Thin-plate spline 


WT 


Wavelet Transform 
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